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Abstract 

The  goal  of  this  project  was  to  use  real  time  data  from  the  US  Naval  Academy's 
Prototype  Communications  Satellite  (PC-Sat)  to  calculate  atmospheric  density  in  the  satellite's 
orbit  over  various  time  intervals.  This  involved  using  the  latitude,  longitude,  altitude,  and  time 
data  from  a  GPS  receiver  on  board  PC- Sat  and  transforming  them  into  the  orbiter's  classical 
orbital  elements  (COEs).  From  these,  the  change  in  the  size  of  the  orbit  can  be  determined  via 
the  change  in  the  semi- major  axis.  Changes  to  the  orbit  are  due  primarily  to  the  non- spherical 
Earth  and  the  atmosphere.  Therefore,  by  accounting  for  the  change  in  semi- major  axis  due  to  the 
non- spherical  Earth,  the  researcher  can  conclude  that  the  remaining  change  is  due  solely  to 
atmospheric  density. 

The  ability  to  determine  atmospheric  density  in  a  specific  orbit  by  knowing  only  the 
position  of  the  satellite  and  few  characteristics  of  the  satellite  itself  will  allow  many  small 
satellites  with  GPS  receivers  to  contribute  to  the  collection  of  data  about  the  upper  atmosphere. 
Being  able  to  measure  the  Earth's  atmospheric  density  with  increased  accuracy  will  then  allow 
satellite  orbit  and  fuel  usage  predictions  to  be  much  more  accurate  and  has  the  potential  to  lower 
the  cost  of  missions. 
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The  Naval  Academy  launched  its  first  midshipman-designed  and  built  spacecraft  on 
September  30, 2001  from  Kodiak,  Alaska.  The  satellite  is  Prototype  Communications  Satellite, 
(PC- Sat),  and  its  main  mission  is  to  “provide  mobile  and  handheld  satellite  digital 
communications  for  amateur  satellite  operators  worldwide  using  the  Automatic-Position- 
Reporting- System  (APRS).”1  The  satellite  also  has  a  Global  Positioning  System  (GPS)  receiver 
on  board  that  transmits  the  satellite's  position  and  time. 

MIDN  John  Young  completed  a  Trident  project  in  2001  in  which  he  developed  an 
algorithm  to  calculate  atmospheric  density  using  real-time  GPS  data.  His  goal  was  to  determine  a 
method  for  calculating  the  density  of  the  atmosphere  at  a  specific  point  at  a  specific  time  in  a 
single  satellite's  orbit.  As  more  satellites  with  GPS  receivers  are  used,  a  more  accurate  model  of 
atmospheric  density  can  thus  be  determined. 

The  atmospheric  density  over  several  periods  of  time  in  the  orbit  of  PC-Sat  has  been 
determined  using  the  GPS  data  received  and  the  algorithm  developed  by  MIDN  Young.  The 
GPS  data  are  received  as  latitude,  longitude,  altitude,  and  time  in  the  Earth- Centered,  Earth- 
Fixed  (ECEF)  rotating  reference  frame  of  the  Earth.  The  satellite  is  most  easily  described  in  a 
reference  frame  which  is  inertial  with  respect  to  the  center  of  the  Earth  and  the  stars.  Therefore 
the  position  in  the  rotating  reference  frame  must  be  transformed  into  the  inertial  reference  frame 
using  rotation  matrices.  Once  the  position  is  known  in  the  inertial  reference  frame,  the  velocity 
vector  of  the  satellite  must  be  calculated  for  each  position  vector.  The  position  and  velocity 
vectors  at  each  time  can  be  used  to  calculate  the  classical  orbital  elements  of  the  satellite  at  each 
time.  The  challenge  to  the  researcher  lies  in  finding  the  velocity  vectors  of  the  satellite  from 
only  the  position  vectors. 
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Until  recently,  GPS  receivers  have  not  been  very  common  on  spacecraft.  Most  users 
have  not  had  a  need  for  the  accuracy  available  with  GPS  for  use  with  space  operations.  The  idea 
of  using  GPS  onboard  spacecraft  is  somewhat  novel,  and  this  is  among  the  first  research  to  use 
GPS  for  a  purpose  other  than  positioning. 


Background 

PC-Sat 

PC- Sat  is  a  10- inch  cube  with  solar  panels  on  all  six  faces.  Five  of  the  six  faces  use  off- 
the-shelf  solar  panels.  The  sixth  face  has  a  solar  panel  designed  for  space  use,  which  failed 
shortly  after  launch.  This  severely  limits  the  amount  of  power  PC- Sat  collects  while  it  is  in  the 
sunlight,  especially  when  the  failed  solar  panel  is  facing  towards  the  sun.  Batteries  must  be  used 
while  the  satellite  is  not  in  view  of  the  sun  and  when  the  satellite  is  using  more  energy  than  it  is 
producing  from  the  solar  arrays.  The  coordinate  system  corresponding  to  PC- Sat  is  shown  in 
Figure  1. 


Figure  1:  PC-Sat  coordinate  system 
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The  orbital  period  of  PC-Sat  is  approximately  100  minutes.  During  its  orbit  it  spends 
a  maximum  of  35  minutes  in  eclipse  and  a  minimum  of  65  minutes  in  the  sun.  Figure  2  shows 
the  eclipse  times  in  minutes  from  October  2001  to  May  2002. 


Figure  2:  PC -Sat  eclipse  times  for  November  2001-May  2002 

The  high  inclination  of  PC- Sat's  orbit  (67  degrees)  and  the  angle  of  Earth's  spin  axis 
sometimes  allow  PC-Sat  to  be  in  the  Sun  for  an  entire  orbit  for  several  days,  the  first  occurrence 
of  which  was  January  10-21, 2002.  During  these  twelve  days  of  continuous  sun,  PC-Sat  had 
enough  power  to  allow  the  GPS  receiver  to  remain  on  continuously.  PC- Sat  does  not  have 
memory  on  board  and  cannot  store  data  to  transmit  once  it  passes  over  a  ground  station. 
Therefore  it  transmits  telemetry  continuously,  but  only  that  which  is  received  by  a  ground  station 
is  recorded.  With  the  help  of  ground  stations  in  South  Africa  and  Great  Britain,  PC-Sat's 
position  according  to  the  GPS  receiver  was  recorded  during  that  period. 


The  data  received  from  the  GPS  receiver  on  PC- Sat  comes  in  the  “GGA”  format 
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shown  in  Figure  3.  It  gives  the  time  the  position  was  taken,  the  latitude,  longitude,  and  altitude 
of  the  satellite,  and  the  number  of  satellites  tracked. 


time  [YRhhmmss.ms]  (2002,  2056.00) 
r~  latitude  [DegMin.Sec]  (37°  47.281 1'  N) 

longitude  [DegMin.Sec]  (144°  58.7417’  W) 


$GPGGA,022056.00 


,  ^747.281  1,n!  1 4459.7417, w!  1,  12,0.7,  803863. 1,M, 0.0, M„*7B 

UJ  I — 


accuracy  of  fix  (l  =  accurate,  0  =  inaccurate) 

number  of  satellites  tracked  (12) _ 

altitude  [m]  (803863.1  m) - 


Figure  3:  GPS  format 


Orbits 

An  orbit  is  an  ellipse,  as  was  discovered  by  Kepler  in  the  1600s.  An  ellipse  is  one  of  the 
curves  in  the  family  of  conic  sections.  Figure  4  shows  the  various  curves  created  by  the 
intersection  of  a  plane  and  a  right  circular  cone,  all  of  which  are  possible  orbits.  The  elliptical 
orbit  is  most  common,  the  circular  orbit  being  a  special  case  of  the  elliptical  orbit.  Parabolic  and 
hyperbolic  trajectories  also  exist,  although  bodies  in  these  orbits  only  pass  by  once. 


Figure  4:  Family  of  conic  sections;  orbital  paths 
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Ellipses  have  foci  as  shown  in  Figure  5.  For  a  satellite  in  orbit  about  Earth,  Earth  is  at  one 
focus  and  the  other  focus  is  vacant.  Perigee  is  the  point  of  the  ellipse  that  comes  closest  to  the 
Earth.  The  apogee  is  the  point  farthest  from  the  Earth. 


perigee 


Figure  5:  Parts  of  an  elliptical  orbit 

A  satellite's  orbit  is  defined  by  six  classical  orbital  elements  (COEs),  much  as  a  position 
on  the  earth  is  described  by  latitude  and  longitude.  The  COEs  are  in  the  inertial  reference  frame 
with  respect  to  the  stars,  meaning  the  coordinate  system  does  not  appear  to  move  relative  to  the 
stars.  In  defining  this  coordinate  system,  the  first  step  is  to  determine  the  first  point  of  aries,  ?, 
which  is  the  line  from  the  Earth  to  the  Sun  on  the  first  day  of  spring.  The  COEs  are  based  on  the 
first  point  of  aries.  Figure  6  illustrates  these  elements. 
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Figure  6:  Classical  Orbital  Elements 

Semi- major  axis,  a,  is  the  length  of  the  long  axis  of  the  ellipse,  and  also  the  mean  radius 
of  the  path.  The  eccentricity,  e,  is  the  variance  of  an  orbit  from  a  perfect  circle,  and  a  measure  of 
the  difference  between  the  long  and  short  axes  of  the  ellipse.  If  the  semi- major  axis  and  semi¬ 
minor  axis,  b,  are  the  same  length,  the  orbit  is  circular  and  the  eccentricity  is  zero.  As  the  semi¬ 
major  axis  increases  relative  to  the  semi- minor  axis,  the  orbit  becomes  more  elliptical  and  the 
eccentricity  increases.  The  eccentricity  of  an  ellipse  is  between  zero  and  one.  When  the 
eccentricity  equals  1,  the  orbit  is  the  special  case  of  a  parabola.  A  hyperbola  has  an  eccentricity 
greater  than  1. 

The  orbit  of  a  satellite  does  not  necessarily  lie  in  the  plane  of  the  Earth's  equator,  and 
inclination,  i,  is  the  angular  measurement  from  the  equatorial  plane  to  the  plane  of  the  orbit.  An 
orbit  is  also  defined  by  the  longitude  at  which  the  plane  of  the  orbit  and  the  equatorial  plane 
cross.  This  actually  occurs  at  two  points,  so  the  crossing  is  defined  as  the  longitude  of  the 
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ascending  node  (O)  and  is  the  longitude  at  which  the  satellite  ascends  from  the  southern 
hemisphere  to  the  northern  hemisphere.  The  angle  from  this  line  to  the  point  of  perigee  is  called 
the  argument  of  perigee  (? ).  The  position  of  the  spacecraft  along  the  orbit  from  the  point  of 
perigee  is  the  true  anomaly  (?)  and  is  measured  as  the  angle  between  the  vectors  connecting  the 
center  of  the  earth  to  perigee  and  the  spacecraft.  These  six  elements,  a,  e,  i,  Cl,  CD,  and  v  define  a 
single  orbit  and  the  satellite's  position  in  that  orbit.  Table  1  shows  the  approximate  values  of  the 
COEs  for  PC- Sat  as  taken  from  NORAD. 


Table  1:  PC-Sat  Classical  Orbital  Elements 


COE 

value 

a 

7178  km 

e 

0.0005 

i 

67  deg 

n 

varies  (-210  deg  in  Jan  2002) 

CD 

varies  (-275  deg  in  Jan  2002) 

V 

varies  with  time 

GPS 

The  Global  Positioning  System  (GPS)  is  made  up  of  3  segments:  the  space  segment, 
control  segment,  and  user  segment.  The  space  segment  is  a  constellation  of  24  satellites  in  semi- 
synchronous  orbit,  meaning  they  each  orbit  the  Earth  twice  a  day.  From  any  point  on  the  surface 
of  the  Earth,  at  least  4  satellites  are  in  view  at  any  time.  This  means  that  a  receiver  on  Earth, 
barring  possible  blockage  by  local  obstructions,  can  determine  its  position  anywhere.  Each  GPS 
satellite  transmits  2  microwave  carrier  signals  with  a  navigation  message  (time  corrections,  space 
vehicle  orbital  data  sets,  etc.)  and  code  signals  (the  fingerprint  of  each  space  vehicle).  The 
control  segment  is  the  ground  network  that  tracks  the  GPS  satellites  and  corrects  the  clocks  on 
each  satellite  to  ensure  very  precise  location  measurements.  The  user  segment  consists  of  the 
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GPS  receivers  and  the  user  community.  The  GPS  receiver  uses  the  time  difference  between 
the  times  sent  from  the  satellite  and  the  receiver  time  to  calculate  what  is  known  as  the  pseudo¬ 
range,  the  distance  from  each  individual  satellite  to  the  GPS  receiver.  A  receiver's  position  is 
computed  using  triangulation  from  four  pseudo-ranges  from  four  different  GPS  satellites.  Figure 
7  shows  the  geometry  of  the  satellites  relative  to  the  receiver,  the  lines  connecting  each  showing 
the  pseudo-range.3  As  more  satellite  signals  are  processed  and  as  the  angles  between  the 
satellites  increase,  the  position  determined  by  the  GPS  receiver  is  more  accurate. 


Figure  7:  Optimum  geometry  of  GPS  satellites  relative  to  receiver  for  position  fix 

For  a  receiver  on  the  ground,  only  satellites  that  are  above  the  horizon  are  in  view  and  are 
available  for  determining  position.  The  actual  number  in  view  is  a  minimum  of  4,  but  could  be  5 
or  6  satellites.  At  800  km,  however,  there  is  an  average  of  thirteen  satellites  in  view  at  any  given 
time.  The  receiver  on  PC-Sat  can  track  only  twelve  of  the  satellites  in  view,  and  on  average 
tracked  ten  to  eleven  satellites  at  any  given  time.4 

The  positions  from  GPS  receivers  are  provided  in  latitude,  longitude,  and  altitude  above 
the  geoid  Earth.  The  Earth  is  not  a  perfect  sphere;  it  is  wider  at  die  equator  than  it  is  between  die 
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poles.  It  also  has  a  slight  egg-shape  to  it.  Therefore,  the  position  on  a  geoid  earth  is  more 
accurate  than  the  position  on  a  spherical  earth.  World  Geodetic  System  1984  (WGS-84)  is  the 
datum  or  mathematical  model  of  the  Earth's  shape  used  by  GPS.5  This  is  also  the  datum  used  for 
maps  used  every  day,  and  thus  makes  the  conversion  from  GPS  position  to  a  map  position  very 
simple. 

GPS  Receiver 

Most  GPS  receivers  are  designed  for  ground  use,  where  the  motion  between  the  receiver 
and  the  GPS  satellites  is  relatively  small.  The  velocities  of  these  terrestrial  GPS  receivers  with 
reference  to  the  Earth  are  also  very  small.  For  a  spacecraft  traveling  at  7-8  km/s,  however,  the 
relative  motion  between  the  GPS  satellites  and  the  GPS  receiver  becomes  an  issue.  “The  DLR 
Orion  GPS  Receiver  [onboard  PC-Sat]  was  originally  designed  for  terrestrial  applications,  but 
has  received  numerous  modifications  to  provide  accurate  and  reliable  tracking  under  highly 
varying  signal  dynamics  encountered  in  space  applications.”6  This  modification  includes  an 
orbit  propagator  in  the  receiver  that  predicts  the  satellite's  position  before  the  GPS  signals  are 
received.  From  this  estimated  position,  the  Doppler  shift  of  the  signals  as  observed  by  the 
receiver  is  calculated.  Knowing  the  expected  shift  in  frequencies  allows  the  receiver  to  track  the 
specific  satellites  in  view  rather  than  searching  through  the  wide  Doppler  frequency  band  for  all 
the  GPS  satellites.7  For  a  GPS  satellite  directly  above  the  receiver  relative  to  the  center  of  the 
Earth,  there  is  no  Doppler  shift.  It  is  well  known  within  the  GPS  user  community  that  using 
satellites  directly  overhead  for  a  fix  have  poor  altitude  accuracies.  The  Orion  receiver  does  not 
discard  signals  from  satellites  directly  overhead,  and  could  have  some  error  in  the  calculation  of 


altitude. 
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The  Inter-agency  GPS  Executive  Board  (IGEB)  published  the  US  Government's 
standard  precision  of  GPS  in  2001,  after  the  removal  of  Selective  Availability.  They  publish  a 
horizontal  accuracy  of  13  meters,  and  a  vertical  accuracy  of  22  meters.8  Since  the  satellite  is  in 
space,  its  receiver  does  not  experience  as  much  refraction  of  the  atmosphere  to  cause  errors,  and 
theoretically  should  have  equal  or  better  accuracy.  The  receiver  has  a  mask  angle  of -15 
degrees,  meaning  it  does  not  look  for  satellites  that  are  below  15  degrees  of  the  horizon  of  the 
antenna.  The  manufacturer  of  the  Orion  GPS  receiver  confirmed  the  accuracy  of  the  receiver.  In 
testing  the  receiver  before  it  was  installed  on  PC- Sat,  the  receiver  had  an  accuracy  of  3.4  meters 
in  position  and  0.8  m/s  in  velocity. 9 

The  antenna  used  on  PC- Sat  to  receive  the  GPS  signals  is  a  small  monopole,  quarter- 
wavelength  antenna  located  on  one  comer  of  PC- Sat.  The  altitude  of  the  orbit  in  which  PC- Sat  is 
located  is  far  below  the  altitude  of  the  GPS  transmitting  satellites,  and  is  therefore  still  within  the 
constellation  of  GPS  satellites.  NASA  has  successfully  developed  a  GPS  receiver,  named 
PiVoT,  for  use  in  highly  elliptical  orbits  and  orbits  that  move  beyond  the  orbits  of  the  GPS 
satellite. 10  This  receiver  has  received  numerous  improvements  for  detecting  very  weak  and 
highly  vaiying  signals.  The  problems  encountered  when  using  GPS  in  this  higher  orbit  are  much 
greater  than  those  experienced  at  the  800  km  altitude  of  PC-Sat.  Although  a  receiver  for  space 
use  does  require  some  modifications,  the  receiver  that  is  on  board  PC- Sat  does  not  need  the  same 
rigorous  changes  that  PiVoT  does. 

Based  on  the  testing  of  the  Orion  receiver,  the  GPS  data  returned  should  have  been  highly 
accurate  and  reliable.  The  tumbling  of  the  satellite,  however,  does  limit  some  reception  to  the 
single  receiving  antenna.  The  antenna  was  placed  on  the  comer  so  that  it  would  have  the  least 
amount  of  area  blocking  it  from  incoming  signals.  Shortly  after  launch  the  GPS  was  turned  on, 
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and  the  receiver  was  unable  to  get  a  single  position.  Several  months  later,  the  receiver  was 
turned  on  again  and  data  were  taken  over  12  days.  The  difference  between  the  two  data- 
gathering  attempts  was  the  tumbling  of  the  satellite.  Between  the  on-board  magnets  aligning  the 
z-axis  with  the  Earth's  magnetic  field  and  the  painted  antennas  regulating  the  speed  of  rotation, 
PC-Sat  settled  into  a  slower  tumble.  This  allowed  the  GPS  receiver  to  pick  up  GPS  signals  and 
predict  forward  to  the  next  expected  signal  frequency.  As  the  satellite  tumbles,  however,  it  still 
could  lose  one  satellite,  pick  up  a  new  satellite,  and  not  realize  what  occurred.  The  receiver 
would  report  an  inaccurate  position,  but  label  it  as  accurate.  Figure  3  shows  how  the  accuracy  of 
the  fix  is  reported.  There  are  several  points  in  the  GPS  data  where  it  appears  this  has  occurred. 
The  positions  of  the  satellite  track  along  an  apparent  orbit,  one  to  twelve  consecutive  positions 
appear  drastically  out  of  place,  and  then  the  positions  return  to  the  expected  orbit.  In  the  12  days 
of  data  and  over  1600  positions  reported,  this  occurs  only  7  apparent  times.  Figure  8  shows  one 
example  of  the  outliers.  The  points  that  are  stars  should  be  between  the  squares  as  expected  for  a 
nearly  circular  orbit 
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NORAD 

"The  North  American  Aerospace  Defense  Command  (NORAD)  is  a  bi-  national  United 
States  and  Canadian  organization  charged  with  the  missions  of  aerospace  warning  and  aerospace 
control  for  North  America.  Aerospace  warning  includes  the  monitoring  of  man-made  objects  in 
space."11  "The  Space  Control  Center  (SCC)  supports  the  US  Space  Command's, 
(USSPACECOM's),  missions  of  surveillance  and  protection  in  space.  The  center’s  primary 
objective  in  performing  the  surveillance  mission  is  to  detect,  track,  identify,  and  catalog  all  man¬ 
made  objects  orbiting  the  earth.  The  SCC  currently  tracks  over  8,300  objects  including 
payloads,  rocket  bodies  and  debris.  Knowing  where  these  objects  are  contributes  to  several 
mission  areas,  including  collision  avoidance  for  the  space  shuttle  crew."12 
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NORAD  is  able  to  track  these  Earth  satellites  using  a  system  called  Space  Detection 
and  Tracking  System  (SPADATS).  This  system  is  "capable  of  detecting  and  tracking  space 
vehicles  from  earth  and  reporting  the  orbital  characteristics  of  these  vehicles  to  a  central  control 
facility."13  SPADATS  is  made  up  of  two  parts:  US  Navy's  Space  Surveillance  System 
(SPASUR)  and  the  US  Air  Force's  SPACETRACK  system.  SPASUR  detects  and  determines  the 
orbital  elements  of  man-made  objects  orbiting  Earth  by  using  an  electronic  fence,  a  continuous 
wave  of  energy  beamed  vertically  across  the  continental  United  States.  This  'fence'  stretches 
from  Georgia  to  California,  1,600  km  off  each  coast,  and  24,100  km  up  into  space.  It  uses 
bistatic  radar  made  up  of  three  transmitters  (each  up  to  two  miles  in  length)  and  six  receivers. 

The  receivers  are  located  separately  from  the  transmitters  and  receive  signals  from  the 
transmitters  that  have  reflected  off  of  an  object  in  space.  Observations  from  two  or  more 
receivers  are  used  to  calculate  the  position  of  the  object. 14  "SPACETRACK  consists  of  a 
worldwide  network  of  radars,  space-probing  cameras,  and  communications.  An  operational 
control  center  with  a  central  data-processing  facility  called  the  Space  Defense  Center,  located  at 
NORAD  headquarters,  serves  to  integrate  the  entire  network  of  both  Navy  and  Air  Force 
information."15  The  result  of  all  these  transmitters,  receivers,  radars,  cameras,  and 
communications  is  a  sophisticated  satellite  tracking  system.  NORAD  then  publishes  the  orbital 
information  of  each  object  on  the  World  Wide  Web  in  the  2-line  element  format.  Figure  9  shows 
a  sample  2- line  element  set  and  what  the  different  numbers  represent. 
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Figure  9:  Sample  NORAD  2 -line  elements 


PC-Sat  Details 

PC- Sat  has  no  propulsion  system  on  board  and  is  free  to  rotate.  It  uses  metal  tape 
measures  of  different  lengths  as  antennas  for  the  different  downlink  and  uplink  frequencies. 

Those  that  are  in  the  x-y  plane  are  painted  white  on  one  side  and  black  on  the  other  to  allow  solar 
pressure  to  create  rotation  about  the  z-axis.  As  photons  from  the  sun  hit  the  antennas  that  are 
painted  black,  they  are  absorbed  by  the  antenna  and  momentum  transfer  occurs.  When  the 
photons  hit  the  white  sides  of  the  antennas,  they  are  reflected,  and  the  antenna  receives  a 
momentum  exchange  of  twice  what  the  black  antenna  receives.  The  greater  torque  on  the  white 
antenna  causes  the  spacecraft  to  spin.  Figure  10  shows  the  rotation  about  the  z-axis  created  by 
pressure  in  the  x-  and  y-axes. 


Figure  10:  Rotation  of  PC-Sat  about  the  z-axis 


Theoretically  the  spacecraft  would  continue  to  speed  up  indefinitely.  However,  PC-Sat  is 
operating  in  the  Earth's  magnetic  field.  As  a  conductor  is  moved  through  a  magnetic  field,  a 
voltage  is  induced  in  the  conductor.  The  current  which  is  produced  creates  a  magnetic  field  in 
the  opposite  direction.  Therefore,  as  PC- Sat  is  spinning  through  the  Earth's  magnetic  field,  an 
equilibrium  spin  rate  is  reached  when  the  magnitude  of  the  force  of  the  induced  magnetic  field 
equals  the  magnitude  of  the  force  of  the  sunlight's  momentum  on  the  antennas.  The  spin  rate  is 
also  moderated  by  the  time  in  eclipse.  When  the  sun  is  blocked  from  the  satellite's  view  by  the 
Earth,  it  no  longer  causes  the  satellite  to  spin,  and  the  latter  slows  down. 

PC-Sat  also  has  temperature  sensors  on  each  solar  panel  face.  These  temperature 
readings  are  sent  back  to  Earth  as  a  part  of  the  telemetry  data  continuously  transmitted  by  PC- 
Sat.  The  USNA  ground  station  collects  these  telemetry  data  every  time  PC-Sat  passes  over 
Annapolis,  Maryland.  CDR  Robert  Bruninga,  USNA  satellite  ground  station  director  and  PC-Sat 
main  operator,  has  monitored  these  temperature  telemetry  data  over  the  past  year  and  has 
determined  the  spin  rate  to  be  about  0.6  ipm  while  orbit  is  in  eclipse  season,  and  0.8  rpm  during 
full  sun  periods. 16 
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The  orientation  of  PC- Sat  with  the  magnetic  field  is  accomplished  by  magnets  placed 
in  PC- Sat.  The  +z  face  is  designed  to  point  towards  the  south  pole,  and  the  -z  face  to  point  to 
the  north  pole.  Therefore,  the  satellite  is  always  aligned  with  the  magnetic  field  of  the  Earth.  If 
PC-Sat  were  in  a  polar  orbit,  meaning  it  had  an  inclination  of  90  degrees  and  went  directly  over 
both  poles,  it  would  have  a  tumble  rate  of  2  complete  tumbles  per  orbit.  This  is  not  easy  to 
visualize,  but  think  of  the  satellite  above  the  North  Pole  with  the  -z  face  towards  the  earth.  Then 
at  the  equator,  the  z-axis  is  aligned  with  the  -z  face  now  pointing  to  the  north  and  the  +z  face 
pointing  towards  the  south.  When  the  satellite  is  over  the  south  pole,  the  +z  face  is  now  facing 
the  earth,  but  it  is  again  in  the  same  orientation  as  it  was  at  the  North  Pole.  Therefore  it 
underwent  an  entire  rotation  about  the  x-  or  y-axis  in  half  of  an  orbit  Figure  1 1  illustrates  this 
motion. 


Figure  11:  Rotation  of  magnetically  aligned  satellite  in  polar  orbit 

PC- Sat  is  not  in  a  polar  orbit,  however,  and  therefore  will  not  completely  roll  over  the  z- 
axis.  Instead  it  will  wobble  along  the  z-axis.  The  motion  that  results  from  the  spin  and  wobble 
is  not  a  simple  addition  of  two  rotations  as  a  spinning  top  is.  Instead  it  is  a  complex  tumble,  the 
modeling  of  which  could  be  a  project  of  its  own. 


Atmospheric  Density 
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The  ideal  gas  law  connects  temperature  T,  pressure  p,  and  density  ?,  of  a  gas  of 
molecular  weight  M,  as  shown  in  equation  (1)  where  R  is  the  gas  constant. 

£  =  —  (1) 
p  M 

The  change  in  pressure  with  height  y  is  given  by  the  hydrostatic  equation  in  equation  (2) 
and  is  accurate  up  to  several  hundred  kilometers. 


dp 

-f  =~pg 
dy 


(2) 


These  two  equations  are  the  basis  for  the  exponential  decrease  of  the  density  of  the 
atmosphere  as  a  function  of  height,  y,  the  result  of  which  is  shown  in  equation  (3)  where  h  is  the 
scale  height,  the  altitude  at  which  the  atmospheric  density  is  1/e  the  density  at  the  surface  of  the 


Earth. 
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P  =  Poe[  h  j  (3) 

In  practice,  it  has  been  determined  that  equation  (3)  loses  accuracy  above  approximately  500 
km. 18  This  is  due  mostly  to  the  mean  free  path  of  the  air  molecules  being  larger  than  the  scale 
height,  and  therefore  the  molecules  have  ballistic  paths  rather  than  colliding  with  one  another. 
Above  this  500  km  threshold  a  much  more  complicated  model  includes  the  effects  of  day  and 
night,  the  Sun's  solar  activity,  the  geomagnetic  planetary  index,  and  a  still  unknown  variability  of 
±10  percent.  Accurate  density  predictions  become  extremely  difficult  as  the  number  of  variables 
increases.  Equation  (3)  also  incurs  some  errors  from  the  way  it  was  calculated.  The  acceleration 
due  to  gravity,  g,  is  also  not  a  constant  as  the  altitude  increases,  and  the  simple  integration  of 
equation  (2)  to  equation  (3)  is  not  accurate  over  large  differences  in  height. 
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There  are  many  published  models  of  atmospheric  density,  broken  into  various 
categories  and  taking  into  account  a  range  of  variables.  Some  cover  specific  conditions  or 
geographic  regions,  such  as  the  International  Tropic  Reference  Atmosphere.  Others  are  general 
models  for  the  whole  Earth,  such  as  the  Jacchia  J70  model.  The  choice  of  density  model  is 
dependent  upon  the  purpose  of  its  use,  and  so  the  selection  of  a  model  can  be  difficult. 1 

The  atmospheric  density  that  is  expected  at  800  km  using  the  exponential  model 
(Equation  3)  is  on  the  order  of  10’ 14  kg/nr?. 


Young's  Algorithm 

In  MIDN  John  Young's  Trident  project,  he  developed  an  algorithm  for  calculating 
atmospheric  density  from  a  satellite's  GPS  data.  He  used  Satellite  Tool  Kit  (STK)  and  ran  two 
models  of  a  satellite's  orbit;  one  without  atmospheric  drag  and  one  with  atmospheric  drag.  He 
used  the  difference  between  the  two  models  to  calculate  the  change  in  the  orbit  due  to  drag.  The 
position  and  velocity  data  taken  from  a  GPS  receiver  were  transformed  into  the  COEs  of  the 
orbit.  The  semi- major  axis,  a,  is  a  constant  throughout  an  undisturbed  orbit,  and  therefore  the 
change  in  a  is  a  good  measure  of  the  change  in  the  orbit.  The  density  can  be  related  to  the 
change  in  the  semi- major  axis  over  time,  as  well  as  to  the  eccentricity  of  the  orbit  and  the  true 
anomaly.  Equation  (4)20  shows  the  equation  for  atmospheric  density  from  COEs. 
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The  mass  (m),  coefficient  of  drag  (Cd),  and  cross-sectional  area  (A)  are  all  properties 
of  the  satellite.  The  velocity  of  the  satellite  relative  to  the  atmosphere,  rather  than  to  the  surface 
of  the  Earth,  is  yei.  The  classical  orbital  elements  in  this  equation  are  n,  e,  ?,  and  a. 

Young  made  several  assumptions  in  his  analysis.  One  was  that  the  spacecraft  had  a 
constant  cross-sectional  area  and  coefficient  of  drag.  Another  was  that  the  atmosphere  rotates 
with  the  surface  of  the  earth.  He  also  assumed  that  the  effects  on  the  orbit  due  to  the  shape  of  the 
earth  and  the  atmosphere  were  independent  of  one  another.  Young  concluded  that  the  largest 
potential  source  of  error  would  be  the  rotation  speed  of  the  atmosphere  relative  to  the  satellite. 

To  solve  for  the  expected  change  in  semi- major  axis  over  time  for  a  given  density, 
Equation  4  can  be  worked  backwards.  For  the  density  of  the  exponential  model  of  1x10  14 
kg/m3,  the  expected  change  in  semi- major  axis  over  time  is  -8.9x1 0"6  m/s.  Table  2  shows  the 
expected  changes  over  various  time  intervals. 

Table  2:  Expected  changes  in  semi-major  axis  for  density  of  10A-14  kg/mA3 


Time  interval 

Change  in  semi- major  axis 

100  minute  orbit 

-0.053  meters  j 

1  day 

-0.768  meters  1 

12  days 

-9.212  meters 

1  year 

-280.2  meters 

Spacecraft  Mass 

PC- Sat  has  no  propulsion  system  on  board  and  therefore  has  a  constant  mass.  The  mass 
of  the  spacecraft  is  12  kg.  Therefore  in  the  equation  for  atmospheric  density,  this  mass  is  also  a 
constant  for  all  orientations  and  positions. 
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Cross-Sectional  Area 

In  order  to  calculate  the  atmospheric  density  from  Young's  algorithm,  the  cross-sectional 
area  of  PC-Sat  is  needed.  Because  it  is  tumbling,  the  area  is  continuously  changing.  Therefore 
an  assumed  average  value  of  the  cross-sectional  area  for  the  entire  orbit  is  needed.  This  average 
value  was  calculated  in  a  way  similar  to  how  the  average  power  output  of  the  solar  panels  was 
calculated.  We  know  that  the  rotation  of  PC-Sat  is  random  and  is  a  tumbling  motion  as  opposed 
to  a  spinning  motion.  Assuming  each  orientation  has  an  equal  chance  to  be  facing  the  direction 
of  motion,  each  orientation  was  weighted  according  to  its  frequency  of  occurrence  on  the 
spacecraft.  A  cube  has  6  faces,  12  edges,  and  8  comers  for  a  total  of  26  ‘‘basic”  orientations. 

The  three  integral  orientations  (face,  edge,  and  comer)  are  not  the  only  possible  orientations  as 
the  satellite  can  be  oriented  at  any  angle  between.  Figure  12  shows  these  three  orientations,  their 
dimensions,  and  their  respective  areas. 


Comer  View 


Figure  12:  Cross-sectional  areas  for  the  3  major  orientations 

The  face  view  occurs  six  times  out  of  26  orientations,  so  the  area  of  100  in2  is  weighted 
6/26.  After  weighting  each  of  these  areas  according  to  their  frequency  of  occurrence,  the 
average  cross-sectional  area  is  140.9  in2  (0.091  m2).  A  sphere  touching  the  comers  of  the  cube 
would  have  a  cross-sectional  area  of  157.0  in2  (0.101  m2). 
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NORAD  includes  in  its  2- line  elements  a  term  called  Bstar  (B*).  This  is  not  the 
ballistic  coefficient,  but  is  related  to  it.  The  ballistic  coefficient  is  a  part  of  the  equation  used  to 
calculate  density  from  the  change  in  semi- major  axis.  Equation  (5)  shows  how  the  two  terms  are 
related. 


RePo 
2  B* 


=  BC 


(5) 


The  atmospheric  density  here  is  the  density  at  the  orbit's  perigee  and  is  assumed  to  be 
constant  for  a  given  altitude.  The  radius  of  the  Earth  is  6378.135  km,  giving  a  simplified 
equation  for  the  ballistic  coefficient.21 


BC  =  - 


1 


kg 


12.7416215*  m2 

Figure  13  shows  the  variance  of  the  ballistic  coefficient  over  a  year.  The  points  that  are 


(6) 


square  are  the  days  for  which  GPS  data  were  collected  in  January  of  2001. 
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Figure  13:  B*  values  from  NORAD  2  line  elements 


The  method  by  which  NORAD  calculates  this  B*  value  is  not  well-published.  In  order  to 
determine  the  feasibility  of  using  a  constant  ballistic  coefficient,  a  comparison  between  the 
atmospheric  densities  calculated  using  the  variable  B*  and  the  constant  ballistic  coefficient  was 
done.  The  result  is  shown  in  Figure  14.  NORAD  apparently  uses  a  constant  density,  based  on 
the  altitude  of  the  satellite,  to  calculate  a  B*  to  fit  a  smooth  atmospheric  model.  Figure  14  does 
not  show  this  smooth  density  that  is  expected  with  the  B*  values,  but  it  does  show  a  smooth 
density  for  the  constant  ballistic  coefficient.  The  scale  of  the  graph  is  also  important  to  note,  as 
the  density  is  of  the  same  magnitude  as  predicted  using  equation  (3).  The  values  for  atmospheric 
density  that  were  calculated  using  equation  (4)  and  NORAD  data  only  are  reasonable. 
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NORAD  density  over  time 
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Figure  14:  Comparison  of  variable  B*  and  constant  ballistic  coefficient 


Aerodynamic  Drag :  Coefficient  of  Drag 

It  is  important  to  know  the  density  of  the  atmosphere  through  which  objects  are  traveling 
because  of  the  aerodynamic  drag  force  it  creates  on  the  object.  Although  spacecraft  operate  in 
space,  there  is  still  a  significant  amount  of  atmosphere  at  low- earth  orbit  (LEO)  altitudes.  At  an 
altitude  of  800  km,  the  atmospheric  density  is  approximately  twelve  to  thirteen  orders  of 
magnitude  less  than  that  found  at  sea  level,  around  10'14  kg/m3.22  This,  however,  is  still 
significantly  denser  than  the  vacuum  of  space  at  higher  altitudes  and  therefore  this  drag  effect 
must  be  taken  into  account  at  lower  altitudes.  The  amount  of  drag  an  object  experiences  is 
dependent  on  the  amount  of  atmosphere  present,  or  the  atmospheric  density.  Therefore,  the  drag 
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force  due  to  the  atmosphere  must  be  accounted  for  in  predicting  the  orbit  of  a  satellite  in  low- 
earth  orbit. 

There  is  a  major  difference  between  airflow  on  the  surface  of  the  earth  and  the  flow  of 
the  atmosphere  at  orbital  altitudes.  The  airflow  experienced  on  Earth  as  wind  and  weather  is  the 
result  of  continuum- flow  aerodynamics,  meaning  the  particles  in  the  atmosphere  continuously 
collide  with  one  another  and  interfere  with  the  path  of  other  particles.  Satellites  generally 
experience  free- molecule  flow  in  which  the  molecules  are  reflected  or  re-emitted  from  the 
spacecraft  and  do  not  interfere  with  approaching  or  incident  molecules.  This  is  a  legitimate 
assumption  for  spacecraft  above  200  km  because  of  the  large  mean  free  path  of  the  air  molecules 
relative  to  the  size  of  the  satellite.23 

Aerodynamic  drag  is  created  by  particles  in  the  atmosphere  colliding  with  a  spacecraft 
and  transferring  part  of  the  particle’s  momentum  to  the  spacecraft.  The  exact  amount  of 
momentum  transfer  is  dependent  on  the  angle  of  the  surface  compared  to  the  fluid  flow. 

Equation  (7)  gives  the  momentum  transfer  for  any  collision  where  pi  is  the  momentum  of  the 
incident  particle  and  pr  is  the  momentum  of  the  reflected  particle.  Equation  (8)  shows  equation 
(7)  in  terms  of  pi,  which  is  usually  a  known  quantity. 24  The  angle  ?  is  the  angle  of  incidence 
measured  with  respect  to  the  normal  of  the  surface. 


4 P=  Pi +  Pr 


l  Pi 


=  pi(l+f(e))  =  mv(  1  +  f(6)) 


(7) 

(8) 


In  a  given  amount  of  time,  At,  the  total  mass  deposited  on  a  surface  is  shown  in  equation  (9). 

m  -  pAvAt  (9) 
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By  substituting  equation  (9)  into  equation  (8)  and  then  using  the  definition  of  a  force  as  the 
change  in  momentum  over  the  change  in  time,  the  result  is  equation  (10). 

Fd  =±pv2„lACD.  (10) 

The  atmospheric  density  is  ?,  the  velocity  of  the  spacecraft  relative  to  the  atmosphere  is  vrei,  the 
cross-sectional  area  is  A,  and  Cd  is  the  coefficient  of  drag  of  the  spacecraft. 

The  coefficient  of  drag  is  a  dimensionless  quantity  used  to  describe  the  interaction 
between  a  fluid  and  an  object  traveling  through  the  fluid.  It  is  defined  as  shown  in  equation 
(11),  taken  from  equations  (8)  and  (10). 

C„  =  2(l+/(0))  (11) 

The  exact  value  of  f(?)  is  difficult  to  calculate  theoretically  because  of  the  uncertainty  of  the 
nature  of  individual  atomic  collisions.  Particles  may  scatter  either  elastically  or  they  may  adhere 
to  the  surface  until  thermal  equilibrium  is  achieved  and  scatter  randomly.  The  first  is  called 
specular  reflection;  the  latter  diffuse  reflection. 25  Therefore  the  most  accurate  method  of 
determining  a  spacecraft's  coefficient  of  drag  is  to  determine  it  experimentally.  Drag  coefficient 
values  can  range  from  1.9  to  2.6,  with  the  average  value  at  2.20. 26 

The  coefficient  of  drag  is  constant  for  a  spacecraft  that  has  a  constant  orientation  with 
respect  to  the  atmosphere.  This  is  true  for  satellites  that  are  perfect  spheres  or  have  attitude 
control  systems  that  keep  them  in  a  specific  orientation.  PC- Sat  does  not  meet  either  of  these 
qualifications,  and  therefore  does  not  have  a  constant  coefficient  of  drag.  It  is  tumbling  through 
space  and  the  faces  of  the  cube  are  constantly  changing  angles  relative  to  the  atmosphere.  Figure 
15  shows  experimental  coefficients  of  drag  for  a  flat  plate  at  various  angles  relative  to  the 
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direction  of  motion. 
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Figure  15:  Experimental  Coefficients  of  Drag  for  a  flat  plate  at  varying  angles 


If  PC- Sat  were  flying  face-on  with  the  normal  of  the  face  in  the  direction  of  motion,  it 
would  have  a  coefficient  of  drag  of  about  2.5.  If  it  was  flying  with  an  edge  in  the  direction  of 
motion,  PC- Sat  would  effectively  be  two  flat  plates  flying  at  a  45  degree  angle  and  the 
coefficient  of  drag  would  be  about  2.00.  If  PC-Sat  had  a  comer  in  the  direction  of  motion,  it 
would  effectively  be  three  flat  plates  flying  with  an  angle  of  incidence  of  45  degrees.  This  would 
also  give  a  coefficient  of  drag  of  about  2.00.  Using  the  same  method  as  for  the  cross- sectional 
area  to  weight  the  different  occurrences  of  each  orientation,  the  average  coefficient  of  PC-Sat  is 
2.11. 

The  error  that  comes  from  using  the  average  value  of  2.1 1  for  the  drag  coefficient  can  be 
estimated.  The  equation  for  density  as  given  in  equation  (4)  shows  that  density  is  inversely 
proportional  to  the  coefficient  of  drag.  If  the  average  coefficient  of  drag  is  2. 1 1  and  the 
instantaneous  coefficient  of  drag  is  the  maximum  theoretical  value  of  2.5,  the  increase  in  drag  is 
1 8  percent  of  the  average.  This  in  turn  causes  a  change  in  density  of  about  1 8  percent. 
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Considering  the  goal  of  using  this  method  is  to  reduce  the  uncertainty  in  measurement  to 
below  15  percent,  using  an  average  value  for  coefficient  of  drag  will  not  allow  this.  The 
conclusion  section  will  further  address  this  issue. 

Relative  Velocity 

The  velocity  of  the  satellite  relative  to  the  fluid  through  which  it  is  traveling  is  the 
relative  velocity  referred  to  in  equation  (4).  The  velocity  that  is  calculated  for  the  satellite  in  the 
inertial  reference  frame  is  relative  to  the  center  of  the  earth.  The  velocity  of  the  atmosphere 
relative  to  center  of  the  center  of  the  earth  must  be  known  in  order  to  calculate  the  velocity  of  the 
satellite  relative  to  the  atmosphere. 

King-Hele  derives  an  equation  for  the  velocity  relative  to  the  atmosphere  dependent  on 
the  spacecraft's  velocity  relative  to  the  center  of  the  Earth,  v,  the  spacecraft's  distance  from  the 
center  of  the  Earth,  r,  and  the  inclination  of  the  spacecraft's  orbit,  i.28  Equation  (12)  shows  the 
result,  with  w  as  the  angular  rotation  of  the  atmosphere  in  the  east- west  direction. 

vre,  =  v|l  -  —cos  i  j  (1 2) 

If  the  atmosphere  is  assumed  to  rotate  with  the  surface  of  the  Earth,  at  800  km  the 
atmosphere  will  be  moving  at  approximately  0.4  km/s.  The  atmospheric  rotation  in  the  north- 
south  direction  is  assumed  to  be  zero  because  it  is  significantly  smaller  than  the  east- west  winds 
of  0.4  km/s.  PC-Sat  has  an  orbital  speed  on  the  order  of  7.5  km/s.  The  east-west  winds  are  then 
seen  to  be  about  five  percent  of  the  orbital  speed  In  equation  (4),  the  atmospheric  density  is 
inversely  proportional  to  the  square  of  this  relative  velocity.  Altering  the  relative  velocity  by 
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five  percent  gives  a  change  in  the  density  of  about  ten  percent.  If  this  were  the  only 
uncertainty  in  the  calculation  of  atmospheric  density,  this  would  give  a  more  accurate  density 
than  the  current  model. 


Velocity  from  Position 

The  GPS  data  received  from  PC-Sat  are  in  the  form  of  time,  latitude,  longitude,  and 
altitude.  The  vital  piece  of  information  that  is  missing  is  the  velocity  vector.  Knowing  the 
position  vectors  of  the  satellite  at  two  different  times  should  make  the  velocity  easy  to  find 
Unfortunately  this  is  not  the  case.  The  path  over  which  the  satellite  has  traveled  during  the  time 
between  GPS  fixes  is  not  a  straight  line,  so  v  =  ?d/?  t  is  not  an  accurate  calculation.  Figure  16 
shows  the  difference  in  distances  between  the  straight  and  curved  sections.  This  difference  leads 
to  large  errors  in  the  calculated  velocity.  The  velocity  of  the  satellite  is  very  important  because  it 
is  necessary  for  extracting  the  classical  orbital  elements  from  the  GPS  positions,  the  COEs  that 
are  in  turn  necessary  to  calculate  the  atmospheric  density. 


Figure  16:  Difference  in  path  length  of  curve  and  straight  line 


GPS  and  NORAD 


33 


One  source  of  PC- Sat’s  velocity  is  the  NORAD  data.  The  NORAD  classical  orbital 
elements  are  taken  only  once  or  twice  a  day.  However,  the  velocity  can  be  easily  determined 
from  the  classical  orbital  elements  for  those  one  or  two  readings  each  day.  Figure  17  shows  the 
velocities  from  the  NORAD  data  through  the  year. 


Figure  17:  Velocity  of  PC-Sat  from  NORAD  data 

Although  the  velocities  appear  to  be  randomly  scattered,  they  are  in  a  very  narrow  window 
between  7451  m/s  and  7456  m/s.  The  reason  the  velocity  is  changing  over  time  is  because  the 
NORAD  elements  are  averages  over  a  section  of  the  orbit.  The  orbit  is  slightly  elliptical,  so  the 
velocity  at  perigee  is  slightly  faster  than  the  velocity  at  apogee.  The  NORAD  elements  are  taken 
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at  the  same  position  from  Earth  each  time,  but  the  satellite  is  in  a  different  part  of  its  orbit, 
and  so  the  velocity  captured  by  NORAD  changes  slightly.  This  average  velocity  is  not  suitable 
for  calculations  of  atmospheric  density.  The  instantaneous  values  are  necessary  for  any  analysis 
of  the  atmospheric  density  from  GPS  data. 

Herrick-Gibbs  Method 

One  method  of  calculating  instantaneous  velocity  is  to  use  only  position  vectors  to 
calculate  the  orbit  and  thus  the  velocity.  This  Gibbs  method  uses  three  time- sequential  coplanar 
position  vectors  of  a  satellite  in  its  orbit.29  The  result  of  the  method  is  a  velocity  vector  that 
corresponds  to  the  middle  position  vector.  When  this  method  is  used  the  velocity  ranges  wildly 
with  an  apparent  maximum  velocity  about  7.6  km/s.  This  indicates  that  the  method  used  to 
calculate  velocity  is  incorrect  rather  than  the  data.  The  lower  velocities  could  be  due  to  larger 
time  spans  between  position  vectors  and  thus  a  similar  result  as  shown  in  Figure  16  above. 

The  Herrick-Gibbs  method  is  similar  to  the  Gibbs  method  in  that  it  uses  3  consecutive 
position  vectors,  but  it  also  incorporates  the  time  between  the  fixes  into  the  solution  for  the 
velocity.  The  method  uses  a  Taylor-series  expansion  to  obtain  the  velocity  vector  for  the  middle 
position  vector.  The  method  is  specifically  designed  for  position  vectors  taken  by  a  ground 
station,  which  are  generally  very  close  together.  This  is  equivalent  to  the  data  received  from  the 
GPS  receiver,  as  they  are  also  position  vectors  spaced  close  together.  Figure  18  shows  the 
velocities  calculated  using  Herrick-Gibbs  method. 

The  Herrick-Gibbs  program  includes  a  test  to  ensure  the  three  data  points  are  actually  co¬ 
planar  before  the  velocity  vector  is  calculated.  Because  PC- Sat  does  not  have  memory  on  board 
to  store  data  between  passes  over  the  ground  station,  only  the  data  retrieved  real-time  are 
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available.  The  result  is  gaps  in  the  data  between  ground  stations  across  which  the  velocity 
vectors  cannot  be  accurately  calculated.  Therefore,  if  the  spacing  between  two  of  the  three 
position  vectors  is  greater  than  0.3  degrees,  the  velocity  vector  is  set  to  zero. 


The  majority  of  the  velocity  magnitudes  are  around  7.45  km/s  as  is  expected  from  the 
NORAD  average  velocities.  The  mean  velocity  calculated  using  the  Herrick  Gibbs  method  is 
7.4484  km/s.  From  Figure  18,  any  velocity  above  7.47  and  below  7.44  appears  to  be  an  outlier. 

If  these  points  are  not  included,  the  mean  velocity  is  7.4528  km/s.  There  appear  to  be  a  large 
number  of  outliers  on  this  plot.  When  examined  closer,  these  points  are  almost  all  the  first  or  last 
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few  GPS  positions  in  a  group.  This  phenomenon  will  be  explained  later  in  the  GPS  Outliers 
section. 


Kalman  Filter 

An  alternate  way  to  calculate  the  instantaneous  velocity  vector  is  using  a  Kalman  filter. 
This  is  “a  technique  for  computing  the  best  estimate  of  the  state  of  a  time-varying  process.”  It 
is  a  predictor-corrector  technique,  where  the  state  is  estimated  at  each  successive  observation 
time  to  compare  to  the  actual  observation.  The  Kalman  filter  also  “carries  all  the  information 
concerning  past  measurements  in  its  current  state  and  covariance  estimates  and  therefore  doesn’t 
need  to  reprocess  all  of  the  past  measurement  information  at  each  step.”31  Because  an  orbit  is 
not  linear,  an  extended  Kalman  filter  must  be  used.  The  extended  filter  has  a  two-step  method. 
The  new  state  estimate  is  predicted  using  previous  data,  and  then  the  prediction  is  updated  with 
the  new  observations.  It  uses  a  set  of  matrices  to  bias  the  calculation  of  the  propagated  set  of 
data  based  on  accuracies  of  the  current  set  of  data. 

The  extended  Kalman  filter  requires  a  starting  position  and  velocity  vector.  The  NORAD 
data  are  very  useful  for  this.  A  Matlab  program  titled  "NORADtoRVlO.m"  transforms  the  COEs 
from  the  2-line  elements  of  January  10,  2002  into  position  and  velocity  vectors  measured  at  the 
time  of  the  NORAD  reading.  This  is  the  first  NORAD  data  point  after  the  start  of  the  GPS  data. 
The  Matlab  program  "Predict.m"  turns  the  position  and  velocity  vectors  back  into  the  COEs  and 
propagates  them  forward  to  the  next  observation  time.  The  predict  program  uses  only  the  effect 
due  to  the  oblateness  of  the  Earth  and  atmospheric  drag  to  predict  the  future  orbital  elements. 

The  predicted  position  is  compared  to  the  GPS  measured  position  and  the  difference  is  recorded 
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in  a  matrix.  The  bias  matrices  then  are  used  to  calculate  a  final  and  "filtered"  position  and 
velocity  vector  using  both  the  predicted  and  measured  data.  Figure  19  and  Figure  20  show  the 
velocities  calculated  using  the  Kalman  filter. 


Figure  19:  Velocity  from  Kalman  filter:  same  scale  as  HGibbs  velocity 
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Figure  20:  Velocity  from  Kalman  filter:  different  scale  from  HGibbs  velocity 


Orbit  Prediction:  Non- Spherical  Earth 

The  prediction  program  used  for  the  Kalman  filter  only  takes  two  disturbances  into 
account  as  it  propagates  the  satellite’s  position  forward  The  orbit  of  a  satellite  would  be  a 
perfect  ellipse  if  it  occurred  about  a  point  mass  in  a  perfect  vacuum.  The  classical  orbital 
elements  would  be  constants  and  the  satellite  would  return  to  the  same  position  in  space  every 
oibit.  The  earth  is  not  a  point  mass,  and  space  is  not  a  perfect  vacuum,  especially  at  800  km 
where  PC- Sat  is  orbiting.  The  different  perturbations  affecting  orbits  can  be  accounted  for 
mathematically. 

"The  Earth  is  appreciably  oblate,  the  equatorial  diameter  being  42.77  km  greater  than  the 
polar  diameter.'32  The  result  of  this  oblateness  is  a  change  in  the  right  ascension  of  the 
ascending  node,  Q,  and  a  change  in  the  argument  of  perigee,  0).  The  orbit  rotates  about  the  Earth 
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in  the  inertial  reference  frame,  and  the  orbit  rotates  in  the  orbital  plate.  The  amount  of  each 
rotation  is  each  dependent  on  the  semi- major  axis,  the  eccentricity,  and  the  inclination  of  the 
orbit  as  shown  in  equations  (13)  and  (14). 33 

o  (  V 

Ci  =  --nJ2  - 7-  cos i  (13) 

2  2{a(l-e2)J 

,  (  D  f 5 

t b  =  -nJ2  - j-  (5cos2i-l)  (14) 

4  2^(l-e2)J 

Each  of  these  equations  is  a  simplified  version  using  only  the  J2  term,  the  expanded 
version  of  each  uses  Jn  values  2  and  above.  The  numerical  values  for  the  even  Jn  terms  were 
developed  using  accurate  measurements  of  Cl  for  many  satellites  in  different  inclinations  and 
solving  for  the  smaller  Jn  values.  The  odd  Jn  values  were  determined  by  the  same  method  but 
using  the  equation  for  the  periodic  variation  in  perigee  distance.  Table  3  shows  the  numerical 
values  of  J„.34 


Table  3:  Numerical  values  of  Jn 


10yJ2  =  1082626±1 

10yJ5=-245±5 
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10yJ4  =  -1624±2 

10yJ7  =  -336±6 

109Jio  =  -242±9 

Equations  (13)  and  (14)  are  shown  in  the  simplified  form  using  only  the  J2  term  and 
assuming  the  other  Jn  terms  are  sufficiently  small  to  be  disregarded.  The  critical  inclination  for 
the  change  in  longitude  of  the  ascending  node  (equation  (13))  is  90  degrees.  At  this  inclination 
the  longitude  of  the  ascending  node  does  not  change  as  the  cos(90  deg)  is  zero.  Therefore  for 
higher  inclination  orbits  the  change  is  less  than  for  equatorial  orbits.  The  critical  inclination  in 
equation  (14)  is  63.4  degrees,  the  inclination  at  which  there  is  no  rotation  of  the  argument  of 
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perigee  due  to  the  gravitational  field  of  the  Earth.  The  closer  a  satellite's  orbit  is  to  this 
inclination,  the  less  the  Earth's  oblateness  affects  the  position  of  the  argument  of  perigee. 

Orbit  Prediction:  Effects  of  the  Sun  and  Moon 

The  Sun  and  the  Moon  also  affect  orbits  about  the  Earth  by  their  gravitational  attractions 
and  also  by  solar  radiation  pressure  from  the  Sun.  "The  effects  of  the  luni-solar  gravitational 
attractions  are  generally  small  and  periodic:  the  change  in  perigee  distance  rarely  exceeds  2  km 
for  a  close  satellite  (less  than  1500  km)  with  e  <  0.2,  although  much  larger  changes  occur  for 
highly  eccentric  orbits.  The  four  orbital  elements  e,  i,  O,  and  ?  are  all  affected,  but  a  remains 
constant."35  Since  the  focus  of  this  project  is  on  the  semi- major  axis,  a,  the  gravitational  effects 
of  the  Sun  and  Moon  can  be  ignored. 

The  solar  radiation  pressure  is  a  periodic  perturbation,  unlike  the  atmospheric  drag.  As 
referenced  earlier,  this  solar  radiation  pressure  is  the  force  which  causes  spin  about  the  z-axis  on 
PC-Sat,  and  is  not  accounted  for  beyond  this. 

Classical  Orbital  Elements 

The  equation  for  atmospheric  density  (equation  (4))  requires  several  Classical  Orbital 
Elements  (COEs),  including  the  eccentricity,  mean  motion,  true  anomaly,  and  semi- major  axis. 
The  average  values  of  these  COEs  for  PC-Sat  can  be  directly  taken  from  the  NORAD  two-line 
element  sets.  They  can  also  be  calculated  from  the  position  and  velocity  vectors  at  each  instant 


for  which  we  have  data. 


NORAD  Data 
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The  data  from  NORAD  are  initially  in  the  form  of  the  2- line  elements.  Rather  than  the 
specific  position  of  the  satellite,  the  NORAD  data  give  information  about  the  overall  orbit  in 
which  the  satellite  travels.  The  semi- major  axis  is  not  specifically  given  in  the  2- line  element 
set,  but  the  mean  motion  is.  The  relationship  between  mean  motion  (n)  and  semi- major  axis  is 

(15) 

Figure  21  shows  the  semi- major  axis  for  PC- Sat  calculated  from  NORAD  over  10 
months,  with  January  1, 2002  as  day  zero.  The  data  overlaid  are  for  the  days  that  have 
corresponding  GPS  data.  The  decline  in  semi- major  axis  determined  from  the  NORAD  data 
across  the  12  days  of  January  containing  GPS  data  is  60  m.  Compared  to  Table  2,  this  change  in 
semi- major  axis  is  almost  10  times  greater,  which  means  a  factor  often  difference  between  the 
NORAD  density  and  the  exponential  model  density  should  be  expected. 
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Figure  21:  Semi-major  axis  from  NORAD  data 


GPS  Data 

Eccentricity 

The  eccentricity  vector  can  be  determined  directly  from  the  position  and  velocity  vectors 
at  a  point  in  time.  The  magnitude  of  the  vector  is  the  eccentricity  of  the  orbit.  The  direction  of 
the  vector  points  toward  perigee.  PC- Sat’s  orbit  is  almost  circular,  so  the  actual  position  of 
perigee  is  somewhat  difficult  to  define.  Only  the  magnitude  is  necessary  for  final  calculations, 
which  will  not  count  on  the  direction,  but  the  vector  is  needed  to  calculate  the  velocity.  Equation 
(16)  shows  the  equation  for  the  eccentricity  vector,  where  |l  is  the  gravitational  constant. 

v2  V-(r  v)v 
e  - —*■ - - - 


(16) 
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Mean  Motion 

The  mean  motion  is  a  function  of  the  semi-major  axis.  Equation  (17)  shows  the 
relationship. 


n  = 


(17) 


True  Anomaly 

The  true  anomaly  is  the  position  of  the  spacecraft  in  its  orbit,  referenced  to  the  point  of 
perigee.  It  is  calculated  from  the  eccentricity  and  position  vector.  Equation  (18)  gives  the 
equation. 


C0S(v)=liiF 


(18) 


if(r  •  v  <  0),thenv  =  360° — v 


Semi-Major  Axis 

The  semi- major  axis  is  also  a  function  of  the  position  and  velocity  at  any  particular  time. 
Equation  (19)  shows  the  relationship. 


(19) 


Another  method  of  finding  the  semi- major  axis  is  to  use  the  period.  The  period  of  the 
orbit  is  directly  related  to  the  semi- major  axis,  as  shown  in  equation  (20).  The  symbol  p  is  the 
Earth's  gravitational  parameter,  the  product  of  Earth's  gravitational  constant  G  and  Earth's  mass 
m.  The  numerical  value  of  p  is  3.986x1 05  km3/s2. 


P  = 


3/2 


(20) 
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The  difficulty  of  this  method  lies  in  the  fact  that  the  orbit  is  not  a  closed  loop.  The 
satellite  does  not  return  to  the  exact  same  position  every  orbit,  one  of  the  results  of  atmospheric 
drag.  Therefore  a  definition  of  orbital  period  for  PC-Sat  must  be  agreed  upon  before  determining 
the  value  of  the  orbital  period.  The  right  ascension  and  declination  would  be  the  best  coordinates 
system  for  this  because  they  are  inertial  with  respect  to  space  and  the  orbital  plane,  similar  to 
longitude  and  latitude  relative  to  the  surface  of  the  Earth.  Therefore  a  certain  right  ascension 
could  be  defined  as  the  starting  point  of  the  orbit  and  a  period  is  complete  when  the  satellite 
returns  to  that  right  ascension.  The  errors  involved  with  this  calculation  come  from  the  time 
between  data  points.  In  thirty  seconds  PC- Sat  travels  approximately  225  km.  The  time  of 
crossing  can  be  approximated  by  using  linear  interpolation  between  the  two  points  on  either  side 
of  the  reference  right  ascension.  This  method  assumes  that  the  satellite  is  traveling  at  a  constant 
speed  between  the  two  data  points  used  for  interpolation.  This  is  true  for  a  perfectly  circular 
orbit,  but  PC-Sat’s  is  slightly  elliptical. 

Experimentally,  this  method  did  not  give  the  desired  accuracies.  The  reliability  of  the 
semi- major  axis  calculated  greatly  depended  upon  the  right  ascension  or  declination  chosen.  The 
right  ascension  I  chose  was  the  longitude  of  the  Naval  Academy  because  most  of  the  GPS  data 
recorded  were  taken  here,  so  the  position  coordinates  of  PC- Sat  are  much  more  frequent  across 
this  right  declination.  Table  4  shows  the  varying  average  semi- major  axis  as  a  function  of  the 
reference  latitude.  The  difference  between  the  semi- major  axis  calculated  at  the  equator  and  at 
the  top  of  the  orbit  is  about  60  km.  The  chosen  reference  point  for  the  beginning  and  end  of  the 
orbit  has  too  much  influence  on  the  accuracy  to  be  used  with  this  data  set  to  calculate  the 
atmospheric  density  to  within  15  percent.  With  smaller  intervals  between  GPS  positions,  this 
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method  could  be  possible.  The  thirty  second  intervals,  however,  show  in  Table  4  to  be  too 
large  for  accurate  calculation  of  semi- major  axis. 


Table  4:  Semi-major  axis  calculated  using  period 


Reference 
latitude  [deg] 

Mean  semi- major 
axis  [km] 

0 

7134.2 

10 

7079.1 

|  20 

7001.9 

30 

7002.1 

40 

6982.8 

50 

6974.0 

60 

6973.6 

Outliers 

NORAD  Outlier 

The  outlier  in  the  time-series  data  of  the  NORAD  semi- major  axis  (Figure  21)  is 
somewhat  confusing.  We  know  from  orbital  dynamics  that  it  is  physically  impossible  for  the 
satellite  to  have  a  semi- major  axis  decrease  or  increase  that  amount  without  a  propulsion  system 
or  a  significant  collision  with  another  object  in  space.  For  it  to  have  decreased  and  then 
increased  back  to  the  original  orbit  by  two  collisions  is  impossible.  Even  if  it  was  knocked  into  a 
lower  orbit  somehow,  it  would  then  be  traveling  at  a  higher  velocity  and  it  would  not  return  to 
the  same  orbit  it  left.  Therefore  the  single  outlier  in  the  NORAD  data  was  discarded  as  being  a 
corrupted  data  point. 

GPS  Outliers 

The  outliers  in  the  GPS  data  tend  to  occur  infrequently,  several  at  a  time,  and  afterwards 
the  orbit  returns  to  the  expected  orbit.  As  with  the  NORAD  outlier,  the  satellite  cannot  move 
from  one  orbit  to  another  and  back  again  without  some  propulsion  system  or  collision.  PC- Sat 
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does  not  have  a  propulsion  system,  which  leaves  only  collisions  as  a  possibility.  The  return 
of  the  satellite  to  the  predicted  position  after  several  outliers  indicates  that  the  orbit  did  not 
actually  change  as  the  data  suggest. 

The  Kalman  filter  diverged  several  times  rather  than  converging  as  it  was  designed.  The 
divergence  occurred  when  the  predicted  position  and  the  observed  position  were  drastically 
different  than  one  another  and  the  program  could  not  resolve  which  position  was  correct.  The 
source  of  the  divergence  was  GPS  data  points  that  were  not  a  part  of  the  expected  orbit.  Figure 
22  shows  just  one  of  these  outlier  groups  in  the  data,  along  with  the  points  around  them  to  show 
where  the  orbit  is  expected  to  be. 


Figure  22:  Outlier  in  GPS  position  vectors 
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In  Figure  22,  the  small  dots  are  the  satellite's  orbit.  The  larger  squares  clumped 
together  are  the  points  to  either  side  of  the  outliers  and  are  along  the  path  expected  for  the 
satellite’s  orbit.  The  stars  are  the  outliers,  and  should  lie  between  the  two  groups  of  diamonds. 
They  are  in  noticeably  different  positions  from  the  rest  of  the  expected  orbit.  The  prediction 
program  in  the  Kalman  filter  predicted  a  position  vector  to  be  along  the  expected  orbit,  and  then 
compared  the  prediction  to  the  outlier.  The  result  was  a  large  difference  between  the  expected 
and  reported  position  and  filter  was  unable  to  predict  the  next  position,  causing  a  divergence. 

These  few  outliers  are  most  likely  caused  by  the  tumbling  of  the  satellite.  The  GPS 
receiver  that  is  on  PC- Sat  has  in  it  a  program  to  predict  where  each  of  the  GPS  satellites  will  be 
from  one  observation  to  the  next.  Of  the  12  possible  satellites  from  which  the  receiver  can 
receive,  10  to  11  are  constantly  in  view.  The  receiver  tracks  these  satellites  that  are  in  view  and 
predicts  where  they  will  be  for  the  next  pass.  The  tumbling,  however,  is  not  accounted  for  in  the 
program.  The  receiver  rapidly  detects  new  GPS  satellites  in  view,  but  in  these  three  cases, 
probably  is  not  as  quick.  New  satellites  were  most  likely  added  to  the  solution  during  these 
observations,  and  the  receiver  was  either  unable  to  track  them  or  unaware  that  they  were  not  the 
ones  being  tracked.  The  receiver  notes  these  observations  as  being  accurate,  but  the  figures 
above  show  that  they  are  not.  The  Kalman  filter  also  shows  that  when  the  outliers  are  removed 
from  the  data  set,  the  filter  converges  and  produces  positions  and  velocities  at  each  observation 
time.  With  the  outliers  included,  the  Kalman  filter  diverges  at  the  first  one  and  the  program 
stops. 

Sunny  Leung,  German  Space  Operations  Center,  German  Aerospace  Center  (DLR),  has 
also  noticed  these  few  outliers  and  attributes  them  to  the  tumbling  motion  of  PC- Sat.  “Also  bear 
in  mind,  there  is  a  small  spin  on  PC-Sat  which  caused  frequent  change  in  GPS  visibility,  and  this 
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may  lead  to  erroneous  measurements.”36  The  first  and  initial  attempt  to  use  the  GPS  receiver 
was  immediately  after  launch  of  the  spacecraft.  “Rapidly  varying  antenna  attitude  with  respect 
to  the  GPS  constellation  introduced  great  difficulty  to  the  receiver  in  achieving  bit  and  frame 
synchronization  with  a  particular  PRN.  Although  the  Doppler  aiding  algorithm  was  providing 
correct  Doppler  shift  information  to  the  frequency  search  routing  inside  the  receiver,  the 
unfavorable  tracking  condition  caused  difficulty. ...  In  the  second  attempt,  on  the  31  October 
2001,  the  GPS  receiver  was  able  to  track  11-12  satellites. ...  This  observation  indicated  the 
improved  GPS  antenna  attitude  over  the  30  day  period,  which  would  suggest  a  reduced  spin  rate 
and  tumbling  motion  caused  by  the  magnetic  stabilization  onboard.”37  The  tumbling  initially 
caused  massive  problems  for  the  GPS  receiver,  but  the  magnets  in  PC- Sat  have  slowed  the 
tumbling  and  allowed  the  satellite  to  receive  signals  and  calculate  its  position.  However,  the  few 
outliers  that  did  occur  are  most  likely  due  to  the  tumbling  motion  and  can  be  eliminated  from  the 
data  set. 

Results 

Table  5  summarizes  the  average  atmospheric  density  over  the  twelve  day  period  during 
which  the  GPS  receiver  was  activated. 

Table  5:  Summary  of  average  atmospheric  densities  using  various  methods 


Method 

Density  [kg/m5] 

Theoretical 

~lxl0‘l4 

NORAD 

0.91xl0'u 

Herrick  Gibbs 

1.46x1  O' 12 

Kalman  Filter 

1.29x10" 12 

NORAD  Results 
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From  the  NORAD  data,  there  is  a  decrease  in  semi- major  axis  between  each  data  point 
taken,  each  a  half  day  to  a  day  apart.  This  visible  decrease  in  the  orbit  over  the  short  period  of 
time  shows  that  the  atmosphere  at  800  km  is  affecting  the  satellite.  Using  equation  (4)  to 
calculate  the  density  gives  an  average  density  of  0.91x1  O'13  kg/m2  over  the  12  day  period.  The 
average  density  of  the  atmosphere  through  which  the  satellite  traveled  between  NORAD  data 
points  throughout  the  year  2002  was  6.29xl0'14. 

The  predicted  atmospheric  density  between  each  pair  of  NORAD  data  points  can  also  be 
calculated.  Figure  23  shows  the  result  of  this  calculation.  The  density  does  vary  over  each  time 
period,  which  is  expected  because  the  slope  of  the  semi- major  axis  in  Figure  21  is  not  constant. 
The  greater  the  slope,  the  higher  the  atmospheric  density  is  expected  to  be.  The  days  over  which 
GPS  data  were  taken  have  some  of  the  steeper  slopes  in  the  data,  and  are  therefore  expected  to 
have  higher  densities.  There  is  no  known  relationship  between  the  higher  density  recorded 
during  the  time  the  GPS  receiver  was  turned  on  and  the  fact  that  the  receiver  was  on. 
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Figure  23:  NORAD  atmospheric  density  using  constant  BC 


GPS  Results 

The  goal  of  this  project  was  to  be  able  to  calculate  the  atmospheric  density  from  the  GPS 
data  over  near-real  time  intervals.  The  first  calculation  was  the  overall  change  over  the  twelve 
days.  The  change  here  is  expected  to  be  comparable  with  the  change  expected  from  the  NORAD 
two- line  elements,  giving  similar  density  calculations. 

First  used  was  the  Herrick-Gibbs  approach  to  calculate  the  velocity  and  atmospheric 
density.  Once  the  outliers  were  eliminated  (discussed  with  Figure  18),  mean  density  for  the 
twelve  days  was  1.4554x1 0‘ 12  kg/m3.  This  differs  by  an  order  of  magnitude  from  the  simple 
theoretical  density,  as  well  as  the  density  calculated  using  NORAD  data.  The  density  was  also 
calculated  by  using  the  change  in  semi- major  axis  between  each  GPS  data  point,  but  that  method 
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was  not  successful.  The  data  were  not  smoothed  in  any  way,  so  the  calculated  densities  took 
into  account  all  the  noise  in  the  data.  The  result  was  a  mixture  of  positive  and  negative  densities, 
physically  impossible  because  the  density  cannot  be  negative.  The  Herrick-Gibbs  was  not  the 
best  choice  for  this  research.  Had  the  data  points  received  from  the  GPS  receiver  on  PC- Sat  been 
completely  continuous  without  any  large  gaps  in  time,  die  Herrick-Gibbs  method  could  have 
proved  very  useful.  However,  the  number  of  discontinuities  severely  limited  the  applicability  of 
this  method. 

A  better  method  turned  out  to  be  the  Kalman  filter.  Because  it  smoothes  the  data  as  it 
runs,  the  Kalman  filter  gives  results  more  free  of  high  frequency  noise.  It  is  also  able  to  use  all 
the  data  points  in  the  set  without  relying  on  the  points  to  either  side  being  close  in  time.  Using 
the  Kalman  filter,  the  average  atmospheric  density  for  the  twelve  days  was  calculated  to  be 
1.29 19x1  O' 12  kg/m3.  This  is  also  about  two  orders  of  magnitude  different  than  the  theoretical 
and  one  order  of  magnitude  from  the  NORAD  calculated  densities.  The  next  step  was  to 
calculate  the  atmospheric  density  by  binning  the  data  over  vaiying  numbers  of  points.  The 
results  were  also  heavily  dependent  on  the  number  of  points  averaged  into  a  single  calculated 
density.  Table  6  shows  the  mean  of  the  densities  calculated  using  various  numbers  of  data  points 
per  bin. 

Table  6:  Densities  using  Kalman  filter,  averaged  across  various  data  ranges 


Number  of  data 
points  per  bin 

Density 

(kg/m3) 

19 

1.63x10 

219 

1.58xl0‘l2 

519 

1.18x10'° 

819 

7.26x10*° 

1119 

1.19x10*° 

1519 

1.28x1  O' 12 
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The  densities  calculated  using  the  GPS  data  do  not  match  what  they  were  expected  to 
be  using  simplistic  theory,  which  is  not  of  much  concern.  The  theory  uses  an  exponential 
atmospheric  model,  but  it  well  known  that  the  atmosphere  does  not  decrease  exponentially  at  800 
km  altitude.  TheGPS  densities  are  only  about  an  order  of  magnitude  different  from  what  was 
calculated  using  the  NORAD  data.  To  identify  the  source  of  this  difference,  the  data  themselves 
were  examined.  For  an  800  km  orbit  with  an  approximate  eccentricity  of  0.0005,  a  difference  of 
7.2  km  between  perigee  and  apogee  is  expected.  The  plot  of  the  magnitudes  of  the  position 
vectors  is  shown  in  Figure  24.  The  difference  between  the  largest  and  smallest  of  these  numbers 
is  6.9  km. 


Figure  24:  Magnitude  of  position  vectors  from  GPS  data 

The  calculation  of  these  position  vectors  includes  the  eccentricity  of  the  Earth  itself  and 
the  altitude  of  the  satellite.  Without  including  the  eccentricity  of  the  Earth,  the  range  would  be 
19.5  km  rather  than  just  6.9  km.  The  semi- major  axis  calculated  from  these  was  expected  to  be 
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slightly  decreasing  over  time.  Instead  what  they  showed  was  that  the  semi- major  axis  moved 
similarly  to  the  orbits.  Figure  25  shows  that  as  the  satellite  moved  closer  to  the  Earth,  the  semi¬ 
major  axis  seemed  to  decrease,  and  as  the  satellite  moved  away  from  the  Earth,  the  semi- major 
axis  seemed  to  increase. 


Figure  25:  Comparison  of  semi-major  axis  and  position  vector 

The  semi- major  axis  should  not  depend  on  the  satellite's  position  in  the  orbit  and  should 
be  independent  of  the  actual  distance  of  the  satellite  from  the  center  of  the  Earth.  As  the  satellite 
gets  closer,  the  speed  should  increase  and  the  semi- major  axis,  a  function  of  position  and 
velocity  shown  in  equation  (19),  should  remain  constant.  This  is  not  the  case  shown  in  Figure  25 
and  it  appears  that  the  calculated  velocity  is  the  probable  source  of  the  error. 
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I  was  able  to  calculate  a  density  using  the  GPS  data.  The  result  was  different  than 
predicted  by  a  simple  model,  but  this  acceptable  because  the  simple  exponential  model  is  not 
recognized  to  be  correct  at  800  km.  My  results  are  also  only  an  order  of  magnitude  different 
than  the  average  values  calculated  using  NORAD. 


Conclusion  and  Future  Work 

This  project  has  turned  out  to  be  true  research.  The  method  that  was  originally  designed 
to  determine  semi- major  axis  did  not  work,  and  neither  did  several  others.  The  final  result  was 
that  neither  the  Herrick- Gibbs  method  nor  the  Kalman  filter  were  able  to  calculate  the  accurate 
velocities  needed  for  determination  of  the  semi- major  axis.  This  in  turn  did  not  allow  accurate 
calculation  of  the  atmospheric  density. 

A  part  of  the  difficulty  was  that  PC- Sat  was  not  the  optimal  satellite  to  use  for  this 
research.  An  ideal  satellite  would  be  spherical  and,  since  the  effects  of  atmospheric  drag  are  also 
much  greater  at  lower  altitudes,  a  lower  altitude  orbit  would  also  be  better  for  this  analysis.  The 
GPS  receiver  on  PC-Sat  took  position  fixes  about  three  to  four  times  per  rotation,  during  which 
the  coefficient  of  drag  changed.  This  causes  a  significant  error  of  up  to  18  percent.  Another 
difficulty  was  the  lack  of  velocity  data  from  the  GPS  receiver.  The  receiver  can  calculate  these 
data,  but  velocity  was  not  a  part  of  the  data  received  from  PC- Sat.  A  higher  rate  of  positioning 
would  also  give  better  results  in  calculating  density.  The  receiver  on  PC-Sat  was  set  to  fix  at  30 
second  intervals,  during  which  the  satellite  travels  almost  225  km.  Shorter  fix  intervals  would 
give  smaller  errors  between  predictions  and  GPS  positions  and  perhaps  improve  the  accuracy  of 
the  Kalman  filter  in  calculating  the  velocity. 
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Future  work  with  this  particular  data  set  could  continue.  There  are  other  ways  of 
calculating  the  velocity  of  a  spacecraft  given  its  position  than  were  examined  during  this 
research.  It  would  also  be  interesting  to  study  a  different  satellite  with  different  orbital 
characteristics.  A  better  GPS  receiver,  previously  space  tested  and  designed  to  measure  height 
above  the  surface  of  the  Earth,  would  also  improve  the  results  of  this  method  for  calculating 
atmospheric  density.  Being  able  to  use  real-time  or  near-real  time  data  to  calculate  the 
atmospheric  density  for  a  small  interval  of  an  orbit  would  be  a  great  leap  in  knowledge  of  the 
Earth's  upper  atmosphere.  With  current  accuracies  at  or  greater  than  15  percent,  there  is 
certainly  room  for  improvement.  Several  satellites  of  known  ballistic  coefficients  and  circular 
orbits  would  make  this  research  possible. 
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Appendix  A:  Matlab  Programs 


ReadGPS.m 

%z  is  the  length  of  the  arrays 
%data  is  the  [z,4,12]  array  of  all  the  data  taken 
%latS  is  the  array  of  satellite  latitudes 
%longS  is  the  array  of  satellite  longitudes 
%altS  is  the  array  of  satellite  altitudes 
%timeS  is  the  array  of  time  in  [hhmmss] 

%seconds  is  the  number  of  seconds  at  each  time  data  point 
%minutes  is  the  number  of  minutes  at  each  time  data  point 
%hours  is  the  hour  of  the  time  data  poin  t 
%time  is  the  array  of  time  in  [hh  mm  ss] 

file-diskD'; 
if  file— diskD’ 

datal  0=csvread('D:\Text  files\PCSAT__0201 1 0_GPS.txt’); 
datal  l=csvread('D:\Text  files\PCSAT_0201 1  l_GPS.txt’); 
datal  2=csvread(’D:\Text  files' \PCSAT_0201 12_GPS.txt'); 
datal  3=csvread(’D:\Text  files\PCSAT_0201 13_GPS.txt'); 
data  1 4=csvread('D:\T ext  files' \PCSAT_02 01 14_GPS.txt'); 
datal  5=csvread(’ D:\Text  files\PCSAT_0201 15_GPS.txt'); 
datal  6=csvread('D:\Text  files\PCSAT_0201 1 6_GPS.txt'); 
datal  7=csvread(,D:\T ext  files \PCSAT„0201 17_GPS.txt'); 
datal  8=csvread(’D:\Text  files  YPCSAT  02 01 18_GPS.txt'); 
datal  9=csvread('D:\Text  files  \PCSAT_02 01 19_GPS.txt'); 
data20-csvread('D:\Text  files' \PCSAT  _02 01 20_GPS.txf ); 
data2 1 =csvread('D:\T  ext  files\PCS  AT  _0201 21_GPS.txt’); 
end 


z-274; 

dataGPS=zeros(z,4,12); 
[m,n]=size(datalO); 
dataGPS(l:m9:,l)=datalO; 
[m,n]=size(datal  1); 
dataGPS(l:m,:,2)=datal  1; 
[m,n]=size(datal  2); 
dataGPS(l:m,:,3)=datal2; 
[m,n]=size(datal3); 
dataGPS(l:m,:,4)=datal3; 
[m,n]=size(datal4); 
dataGPS(l  :m,:  ,5)=datal  4; 
[m,n]=size(datal5); 
dataGPS(l  :m,:,6)=datal  5; 
[m,n]=size(datal  6); 
dataGPS(l  :m):,7)=datal  6; 
[m,n]=size(datal  7); 
dataGPS(l:m,:,8)=datal7; 
[m,n]=size(datal  8); 
dataGPS(l:m,:,9)=datal8; 
[m,n]=size(data!9); 
dataGP  S(  1  :m,: ,  1 0)=data  1 9 ; 
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[m,n]=size(data20); 
dataGPS(l:nv,ll)=data20; 
[m,n]=size(data2 1 ); 
dataGPS(l:m,:,12)=data21; 


for  i—  1 : 1 2 

latStemp=fix(dataGPS(:,2>i)V100); 

latS(:J:,i)=(latStemp+(dataGPS(:,2,i)-latStemp.*100)/60).*pi/180;  %[rad]  latitude  of  satellite 

longStemp=fix(dataGPS(:,3,i)7100); 

longS(:,:,i)=(longStemp+(dataQ)S(:,3,i)“longStemp.*l  00)/60).*pi./l  80;  %[rad]  longitude  of  satellite 
altS(:,:,i)=dataGPS(:,4,i)/l  000;  %[km]  altitude  of  satellite 

timeSGPS(:,:,i);=:dataGPS(:>lJi);  %[hhmmss]  time  data  was  taken 

hoursGPS(:,:,i)==fix(timeSGPS(:,:,i)./10000); 

minutesGPS(:,:,i)=fix((timeSGPS(:):,i)-hoursGPS(:,:,i).*l  0000)71 00); 
secondsGPS(:,:,i)=fix((timeSGPS(:,:,i)-hoursGPS(:,:,i).*l  0000- minutesGPS(:,:,i).*  100)); 

timeSGPS(:,l,i)=hoursGPS(:,l,i)  *60*60+minutesGPS(:,l,i)  *60+secondsGPS(:)lsi);  %time  in  seconds 

timeGPS(:,l  ,i)=hoursGPS(:,:,i); 

timeGPS(:,2,i)=minutesGPS(:^i); 

timeGPS(:,3,i)=secondsGPS(:5:5i); 

DateGPS(l ,  1  ,i)=i+9;  %days;  10Jan2002=10  etc 

for  j=l:z 

if  hoursGPS(j,l,i)“0  &  minutes GPS(j,l,i)=0  &  secondsGPS(j,l,i)==0 
DaysGPS(j,l,i)=0; 
else 

DaysGPS0,l,i)=DateGPS(l,l,i)+hoursGPSG,U)/24+minutesGPS(j,l>i)/1440+secondsGPSG,l,i)/86400; 

%days 

end 

end 

end 

f=i; 

for  i=l:12 
for  j=l:z 

if  DaysGPSO,:,i)~=0 
Day(f,l)=DaysGPS(j,l}i); 
lat(f,l)=latSQ,l,i); 
long(f,  1  )=longSQ ,  1  ,i); 
alt(f,l)=altSG,l,i); 
hour(f,l  )=hoursGPS(j,l  ,i); 
minute(f,  1  )=minutesGPS(j  ,1  ,i); 
second(f,  1  )=secondsGPS(j ,  1  ,i); 
sec  onds(f,  1  )=timeS  GP  S  (j ,  1  ,i) ; 
time(f,:)=timeGPS(j,:,i); 
f=f+l; 
end 
end 
end 

for  j=2:g 

dtGPS(j,l)=:(Day(j,l)-Day(j-l,l))*86400; 

end 


LST.m 

%LST.m  calculates  the  local  sidereal  time  from  the  julian  day 

thetaG0=  1.75363  7361;  %rad  for  1  Jan2002 

thetaGGPS=thetaG0+1.002737909350795*2*pi*(Day-l);  %rad 
thetaGPS=thetaGGPS+long;  %rad 
thetaGPS=mod(thetaGPSs2*pi); 


GPSPosition.m 


%Read  data  files  into  array  data(z,4,12) 

%turn  data  from  time,  lat,  long,  and  alt  arrays  into  position  in  IJK  frame 
ReadGPS; 

mu-3.98600441 5e5;  %kmA3/sA2 

%Calculate  site  time  data 
LST; 

%ecc  is  the  eccentricity  of  the  earth 
%radius  is  the  radius  of  the  earth 

%altA  is  the  altitude  of  the  ground  point  above  spherical  earth 
%xR  and  zR  are  the  distances  from  the  center  of  the  earth  along 
%  the  x-y  plane  and  in  the  z-direction,  respectively 
%RGroundPoint  is  the  position  of  the  satellite  ground  point  in  XYZ 
%RGPmag  is  the  magnitude  of  the  RGroundPoint  vector 


%ecc=0; 

ecc^O.  08181 922 1 45 6; 
radius=6378.1363;  %km 

xR=(radius./sqrt(l-eccA2.*(sin(lat)).A2)+alt).*cos(lat);  %km 

zR=(radius.*(l-eccA2)./sqrt(l  -eccA2.*(sin(lat)).A2)+alt).*sin(lat);  %km 

radiusE^sqrt(xR  A2+zR  A2); 

rIJKGPS(:,l)=xR.*cos(thetaGPS); 

rIJKGPS(:,2)=xR.*sin(thetaGPS); 

rIJKGPS(:,3)=zR; 

rIJKGPSmag=sqrt(rIJKGPS(:,l)A2+rIJKGPS(:,2)A2-frIJKGPS(:,3)A2); 


NORADtoRVlO.m 


mu=3 .9860044 1 5e5; 
Rearth=6378.1363; 


%kmA3/sA2 

%km 


%read  in  classical  orbital  elements  at  starting  point  from  line  elements 


il=67.0482; 

bigomegal  =2 1 4.3693; 

el=0. 0004963; 

littleomegal  =276.9578; 

Ml=83.0837; 

nl=14.2861444; 

year=2002; 

day  1  =01 0.94044338; 

ndotl  =0.00000955*2; 

ndotl  =ndotl  Sic2*pi/86400A2; 


%inclination,  degrees 
%right  ascension  of  the  ascending  node,  degrees 
%  eccentricity,  no  units 
%argument  of  perigee,  degrees 

%mean  anomaly,  degrees 
%mean  motion,  revs/day 

%epoch  year,  used  for  thetaGO 
%day  and  fraction  of  day  of  tracking  data 
%first  time  derivative  of  mean  motion,  rev/dayA2 

%rad/sA2 


%start  calculations 

nlrad=nl  *2*pi;  %rad/day 

al=(mu/(nlrad/86400)A2)A(l/3);  %km 

%change  degrees  to  radians 
ilrad=il*pi/180; 

bigomegalrad=bigomegal*pi/l  80; 
littleomegalrad=littleomegal  *pi/l  80; 

M 1  r  ad=M  1  *pi/ 180; 

%calculate  eccentric  anomaly  from  mean  anomaly 
ETrad=- 10000; 

Elrad=Mlrad; 

while  abs(Elrad-ETrad)>0.0001 
ETrad=Elrad; 

E 1  rad=Ml  rad-e  1  *  sin(ET  rad); 
end 
El  rad; 

E1=E1  rad*  180/pi; 

%calculate  nu  from  eccentric  anomaly 
nu  1  rad=acos((cos(E  1  rad  )-e  1  )/(l  -e  1  *cos(E  1  rad))); 
if  El  rad  >  pi 
nu  1  rad=2*pi-nu  1  rad; 
end 

nulrad; 

%calculate  r  and  v  in  PQW  coordinates 
pl=al*(l-elA2); 

rlPQabs=pl/(l+el*cos(nulrad)); 
r  1  PQ=[r  1  PQabs*cos(nul  rad);  r  1  PQabs*sin(nul  rad);0]; 
vlPQ=sqrt(mu/(pl)).*[-sin(nulrad);  (el+cos(nulrad));0]; 
vl  PQabs=sqrt(vl  PQ(1  )A2+vl  PQ(2)A2); 

%calculate  in  IJK  coordinates 
R=zeros(3,3); 

R(l,l)=cos(bigomegal  rad)  *cos(littleomegalrad)-sin(big  omega  1  rad)  *sin(littleomegal  rad)  *cos(il  rad); 
R(l,2)=-cos(bigomegalrad)*sin(littleomegalrad)-sin(bigomegalrad)*cos(littleomegalrad)*cos(ilrad); 


%km 

%km,  [P;Q;W] 

%km/s,  [P;Q;W] 

%km/s 


R(1 ,3)=sin(bigomegal  rad)*sin(il  rad); 

R(2,l)=sin(bigomegalrad)*cos(littleomegalrad)+cos(bigomegalrad)*sin(littleomegalrad)*cos(ilrad); 
R(2,2)— sin(bigomegalrad)*sin(littleomegalrad)+cos(bigomegalrad)*cos(littleomegalrad)*cos(ilrad); 

R(2,3)— cos(bigomegal  rad)*sin(il  rad); 

R(3,  l)=sin(littleomegal  rad)*sin(il  rad); 

R(3,2)=cos(littleomega  1  rad)  *sin(il  rad); 

R(3,3)=cos(ilrad); 

R; 

rIJKl=R*rlPQ;  %km 

vIJKl=R*vlPQ;  %km 

rmagl=sqrt(rIJKl  (l)A2+rIJKl  (2)A2+rIJKl  (3)A2);  %km/s 

vmag  1  =sqrt(vl  JK1  (1 )  A2+vI  JK1  (2)A2+vIJKl  (3)A2);  %km/s 

XI  =[rIJKl (l);rIJKl (2);rIJKl  (3);vIJKl  (1  );vIJKl (2);vIJKl (3)];  %initial  state  vector  forkalman  filter 


Predict.m 


mu=3. 98600441 5e5;  %kmA3/sA2 

Rearth=6378.1363;  %km 

J2=0.0010826269;  %assume  units  as  [  ] 

%J2=0; 

%Tum  position  and  velocity  vectors  into  classical  orbital  elements 

h=cross(rI  JK 1 ,  vIJK  1 );  %kmA2/s 

hmag=sqrt(h(  1  )A2+h(2)A2+h(3)A2);  %kmA2/ s 
K=[0;  0;  1  ]; 

LoN=cross(K,h);  %kmA2/s 

LoNmag=sqrt(LoN(l)A2+LoN(2)A2+LoN(3)A2);  %kmA2/s 

e  l=((vmag  1 A2  -mu/rmag  1  )*rl  JKl-dot(rI  JK1  ,vl  JK1  )* vIJK  1  )/mu; 
e  1  mag=sqrt(e  1  (1  )A2+e  1  (2)A2+e  1  (3)A2); 

Energy=vmagl  A2/2~mu/rmagl ;  %kmA2/sA2 


al=-mu/(2*Energy);  %km 

pl=al*(l -elmagA2);  %km 

mml  =sqrt(mu/al  A3);  %rad/s 

il=acos(h(3)/hmag);  %rad 


bigOl  =acos(LoN(l  )/LoNmag);  %rad 

if  LoN(2)<0 
bigO  1  =2*pi-big01 ; 
end 

little01=acos(dot(LoN,el)/(LoNmag*elmag));  %rad 
if el(3)<0 

littleO  1  =2*pi-little01 ; 
end 

nul=acos(dot(el  ,rIJKl)/(elmag*rmagl));  %rad 
if  dot(rIJKl,vIJKl)<0 
nul=2*pi-nul; 
end 

El=acos((elmag+cos(nul))/(l+elmag*cos(nul)));  %rad 
if  nul  >  pi 
El=2*pi-El; 
end 

M 1  =E  1-e  1  mag*sin(E  1 );  %rad 

%Update  for  perturbations,  no  drag  when  ndotl=0 
a2=al-2*al/(3*mml)*ndotl*dt;  %km 
e2=elmag-2*(l~elmag)/(3*mml)*ndotl*dt;  %[] 
big02=big01  -(3*mml  *RearthA2*J2)/(2*pl  A2)*cos(il)*dt;  %rad 
little02-little01+(3*mml*RearthA2*J2)/(4*plA2)*(4-5*(sin(il)A2))*dt;  %rad 
M2=Ml+mml  *dt+ndotl/2*dtA2;  %rad 

p2=a2*(l-e2A2);  %km 

i2=il;  %rad 

%calculate  eccentric  anomaly  from  mean  anomaly 
M2=mod(M2,  2*pi); 

ET— 10000; 

E2=M2; 


while  abs(E2-ET)>0. 00000001 
ET-E2; 

E2=M2+e2*sin(ET); 

end 

%calculate  nu  from  eccentric  anomaly 
nu2=acos((cos(E2>e2)/(l-e2*cos(E2)));  %rad 
if  E2  >  pi 
nu2=2*pi-nu2; 
end 

mm2 =sqrt(mu/a2  A3); 

%Tum  COEs  into  new  position  and  velocity  vectors 

rPQW2=[p2*cos(nu2)/(l+e2*cos(nu2));  p2*sin(nu2)/(l+e2*cos(nu2));  0];  %km 
vPQW2:=:[-sqrt(mu/p2)*sm(nu2);  sqrt(mu/p2)*(e2+cos(nu2));  0];  %km/s 

R=zeros(3,3); 

R(l,l)=cos(big02)*cos(little02)-sin(big02)5*:sin(little02)*cos(i2); 

R(l,2)=-cos(big02)*sin(little02)-sin(big02)*cos(little02)*cos(i2); 

R(l,3)=:sin(big02)*sin(i2); 

R(2,l)=sin(big02)*cos(little02)+cos(big02)*sin(little02)*cos(i2); 

R(2,2)=-sin(big02)*sin(little02)+cos(big02)*cos(little02)*cos(i2); 

R(2,3)^-cos(big02)st:sin(i2); 

R(3}1  )=sin(little02)*sin(i2); 

R(3,2)=cos(little02)*sin(i2); 

R(3,3)=cos(i2); 

rIJK2=R*rPQW2;  %km 
vIJK2=R*vPQW2;  %km/s 

rmag2=sqrt(rIJK2(l)A2+rIJK2(2)A2+rIJK2(3)A2);  %km 
vmag2=sqrt(vIJK2(l)A2+vIJK2(2)A2+vIJK2(3)A2);  %km/s 


X2=[rIJK2(l);rIJK2(2);rIJK2(3);vIJK2(l);vIJK2(2);vIJK2(3)]; 
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Kalmanl  .m 

%Kalmanl  .m  uses  a  kalman  filter  for  the  GPS  data  converted  to  position  vectors. 
%The  initial  position  and  velocity  vectors  are  taken  from  the  NORAD  TLEs. 

mu=3. 98600441 5e5;  %kmA3/sA2 

Rearth-6378.1363;  %km 

%First  read  in  NORAD  TLEs  and  transform  the  COEs  into  position  and  velocity 
%Retums  the  initial  state  vector  XI  ([rx;ry;rz;vx;vy;vz]) 

NORADtoRVIO; 

rIJKl; 

vIJKl; 

rmagl; 

vmagl; 

%Calculate  the  r  vectors  from  the  GPS  data 
%Retums  the  future  state  vector  z  ([rx;ry;rz]) 

GPSPosition; 

%Initialize  H  vector 
H=[l  0  0  0  0  0 
010000 
0  0  1  0  0  0]; 

%Initialize  P  matrix 
P=[l  0  0  0  0  0 
0  1  0  0  0  0 
00  1  0  0  0 
0  0  0  .0001  0  0 
0  0  0  0  .0001  0 
00  0  0  0.0001]; 

%counter  for  filteredX  and  filteredb 
c=l; 

%for  i— 1 : 12 
for  j=l:g 

if  j >=61 9  &  j<=633 
j=634; 
end 

if  j==723 
j=724; 
end 

if  j>=l  190  &  j<— 1 1 92 
j=l 1 93; 
end 

ifj>=1213  &j<=1215 
j=1216; 
end 

if  j>=l 257  &  j<=1258 
j=1259; 
end 

if  j==1085 
j=1086; 
end 

if  j>=l 250  &  j<=1254 


j=1255; 

end 

if  j>— 1 347  & j<=1354 
j=1355; 
end 

%Determine  the  time  difference  between  the  last  data  point  and  the  next 
day2=Day(j,l); 

dt=(day2-dayl)*86400;  %seconds 

%Predict  the  next  position  vector 

%Find  closest  ndot  from  NORAD  data 

%Input  rIJKl ,  vIJKl ,  rmagl ,  vmagl 

%Retum  X2  ([rx;ry;rz])  (rIJK2,  vIJK2,  rmag2,  vmag2) 

Predict; 


%Calculate  phi  matrix 

phi(  1 , 1  )=  1  +(3  *mu*  dtA2*rI  JK1  (l)A2)/(2*nnagl  A5Hmu*dtA2)/(2*rmagl  A3); 
phi(l,2)=(3*mu*dtA2*rIJKl(l)*rIJKl(2))/(2*rmaglA5); 
phi(l  ,3)=(3*mu*dtA2*rIJKl  (1  )*rIJKl  (3))/(2*rmagl  A5); 
phi(2,l)=phi(l,2); 

phi(2,2)=l+(3*mu*dtA2*rIJKl(2)A2)/(2*rmaglA5)-(mu*dtA2)/(2*rmaglA3); 
phi(2,3)=(3*mu*dtA2*rIJKl(2)*rIJKl  (3))/(2*rmagl  A5); 
phi(3,l)=phi(l,3); 
phi(3,2)=phi(2,3); 

phi(3,3)=l+(3*mu*dtA2*rIJKl(3)A2)/(2*rmaglA5)-(mu*dtA2)/(2*nnaglA3); 

phi(4, 1  )=(3  *mu*dt*rl  JK 1  (1  )A2)/(rmagl  A5>(mu*dt)/(rmagl  A3); 
phi(4,2)=(3*mu*dt*rIJKl(l)*rIJKl(2))/(rmaglA5); 
phi(4,3)=(3*rnu*dt*rIJKl(l)*rIJKl(3))/(rmaglA5); 
phi(5,l)=phi(4,2); 

phi(5,2)=(3*mu*dt*rIJKl(2)A2)/(rmaglA5)-(mu*dt)/(rmaglA3); 

phi(5,3)=(3*mu*dt*rIJKl(2)*rIJKl(3))/(rmaglA5); 

Phi(6,l)=phi(4,3); 

phi(6,2)=phi(5,3); 

phi(6,3)=(3*mu*dt*rIJKl(3)A2)/(rmaglA5>(mu*dt)/(rmaglA3); 


phi(l:3,4:6)=[dt  0  0 
0  dt  0 
0  0  dt]; 


phi(4:6,4:6)=[l  0  0 
0  1  0 
0  01]; 

phi; 

%Calculate  predicted  error  covariance 
%Initialize  Q  matrix 
Q=[.01  0  0  0  0  0 
0  .01  0  0  0  0 
0  0  .01  0  0  0 
0  0  0  .0001  0  0 
0  0  0  0  .0001  0 
0  0  0  0  0  .0001]; 


%Calculate  next  Ptemp 
Ptemp=phi*P*phi,+Q; 


%Input  next  measured  position:  z=[rx,ry,rz] 
z  A=  [rl  JKGP  S  (j ,  1  );rl  JKGP  S  (j  ,2);rl  JKGP  S(j,3)]; 

%Calculate  b  vector:  difference  between  predicted  and  measured  position 
b=zA-H*X2; 


%Calculate  Kalman  gain,  K 
R=[.25  0  0 
0  .25  0 
0  0  .25]; 

tempA=H*Ptemp*H’+R; 
temp  B = in  v(  temp  A) ; 
Kgain=Ptemp*H'*tempB; 

%Update  state  change  estimate 
dx=Kgain*b; 


%Update  Error  covariance  estimate 
P=Ptemp  -Kgain*H*Ptemp; 

%Update  state  estimate 
XfmaKX2+dx; 

%print  filtered  position  and  velocity  with  time 
day2; 

Xfrnal; 

r-sqrt(Xfinal(  1 )  A2+Xfinal(2)A2+Xfinal(3)A2); 
v=sqrt(Xfinal(4)A2+Xfinal(5)A2+Xfinal(6)A2); 
r; 
v; 


%save  data 

filteredX(:,c)=Xfinal; 

filteredb(:,c)=b; 

filteredt(:,c)=day2; 

c=c+l; 

%reassign  initial  values  for  next  iteration 

rIJKl=Xfinal(l:3); 

vl  JK 1  =Xfmal(4 : 6); 

rmagl=r; 

vmagl=v; 

dayl=day2; 

end 

%end 


KalmanltoCOE.m 


count=size(filteredX); 

count=count(2); 

Position=sqrt(filteredX(l,:).A2+filteredX(2,:)  A2+filteredX(3,:)  A2)'; 
Velocity=sqrt(fllteredX(4J:).A2+filteredX(5}:).A2+filteredX(65:).A2); 

%calculate  semi-major  axis 

aK(l,:)-l./((2./Position(l,:)>(Velocity(l):)A2/mu)); 

aK=aK'; 

for  j=2:count 

dadtK(l  ,j)=(aK(j,l)-aK(j-l,l))./dtGPS0sl); 
end 

mmK=sqrt(mu./aK  A3); 

hK=cross(filteredX(l:3>:),filteredX(4:6,:)); 

hKabs-sqrt(hK(l,:)A2+hK(2,:)A2+hK(3,:).A2); 

vcrosshK=cross(filteredX(4:6,:),hK); 

fork=l:3 

eK(k,:)=(vcrosshK(kJ:)./mu>(filteredX(k,:)./Position(l>:)); 

end 

eKmag=sqrt(eK(l,:).A2+eK(2,:).A2+eK(3,:).A2); 
for  j=l:count 
if  hKabs(lJ)~=0 

incK(l  j)=acos(hK(3  j)./(hKabs(l  j))); 
else 

incK(ljH); 

end 

nodevectorK(:j)=cross(K,hK(:  j)); 
end 

nodeK(l,:)=sqrt(nodevectorK(l,:)  A2+nodevectorK(2,:).A2+nodevectorK(3}:).A2); 
for  j=l:count 
if  nodeK(lj)~=0 

bigomegaK(  1  j)=acos(nodevectorK(l  j)/nodeK(l  j)); 

littleomegaK(l  j)=acos(dot(nodevectorK(:  j),eK(:  j))/(nodeK(l  j).*eKmag(lj))); 
nuK(l:j)=acos(dot(eK(:j)JfilteredX(l:3j))./(eKinag(l  J).*Position(l  j))); 
else 

bigomegaK(  1  j  )=0; 
littleomegaK(l  j)=0; 
nuK(lj)=0; 
end 

if  nodevectorK(2  j)<0 
bigomegaK(l,j)=2*pi-bigomegaK(l,j); 
end 

if  eK(3  j)<0 

littleomegaK(l  j)=2*pi-littleomegaK(l  j); 
end 

if  dot(filteredX(l :  3  j  ),filteredX(4: 6  j  ))<0 
nuK(lj)=2*pi-nuK(lj); 
end 
end 

bigomegaK=bigomegaK.*  1 80./pi; 
littleomegaK=littleomegaK.*  1 80./pi; 


%filter  out  the  erroneous  data  points:  SMA  below  7170  km  and  above  7186  km 

aboveA=find(aK>7 1 86); 

below  A=fmd(aK<71 70); 

goodA=find(aK>=71 70  &  aK<=71 86); 

count2=size(goodA); 

count2=count2(l); 

aKfmal=aK(goodA); 

Dayfinal=Day(goodA); 

PositionfinaHPosition(goodA); 

Velocityfmal^Velocity(goodA); 
incKfinal=incK(goodA); 
rnmKfmal=mmK(go  odA) ; 
eKmagfmal=eKmag(goodA); 
nuKfinal=nuK(goodA); 

figure; 

scatter  (Day  ,aK); 
hold; 

scatterCDayfinafaKfinaf'g'); 
xlabel('time  [Day]'); 
ylabel('semi~major  axis  [km]'); 
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COEtoDensity.m 

%COEtoDensity.m  takes  the  COEs  from  Kalman  1  toCOE.m  and  calculates  the  density  using  a  running  average 
%over  a  certain  time  period. 

%run  Kalmanl  .m  and  Kalman  1  toCOE.m  first 

%calculate  the  relative  velocity  of  the  satellite  to  the  atmosphere,  with  initial  assumption 
%that  the  atmosphere  rotates  with  the  earth 
omega-7.2921 15856e-5;  %rad/s 
omegaE=[0,0,  omega]; 

Vrel- Velocity  final. *(l-omega.*Positionfinal./Velocityfinal.*cos(incKfinal)); 

Vrel=Vrel.*10A3;  %m/s 

%calculate  the  ballistic  coeeficient 

mass=12;  %kg 

cD=2.11; 

CSArea-140.9;  %inA2 
CSArea=CSArea*0.0254A2;  %mA2 
BC=mass/(cD*CSArea); 

%calculate  the  average  density 
slope  All-poly  fit(Dayfinal,aKfinal,l^ 

density  All=-slopeAll(l).*1000./86400.*(BC).*(l./((mean(Vrel)).A2)).*(mean(mmK  final). *sqrt(l- 
mean(eKmagfmal)  A2)./sqrt(l+mean(eKmagfinal).A2+2.*mean(eKmagfinal).*cos(mean(nuK)))); 

%calculate  the  change  in  semi  -major  axis 
p=i; 

datapoints-1519;  %must  be  an  odd  number;  is  the  number  of  points  to  be  averaged  across 

middle=(datapoints+l)/2;  %this  is  the  point  for  which  the  average  decrease  in  semi-major  axis  will  be 

calculated 

delta-middle- 1; 

order- 1 ; 

for  j -middle  :(count2- middle) 

slope(j,:)-polyfit(Dayfmal((j-delta):(j+delta)Jl),aKfmal((j-delta):(j+delta)Jl),order); 

da(j,l)-0; 

for  loop  -  1:  order 

da(j,l)  -  slope(j,loop)*(order+l  -loop)*Dayfinal(j,l)A(order-loop)+da(j,l); 
end 

Vmean(l  j)-mean(Vrel(l  ,(j-delta):(j+delta))); 
mmKmean(j ,  1  )-mean(mmKfmal((j  -delta) :  (j +delta) ,  1 )) ; 
eKmagmean(l,j)=mean(eKmagfmal(l,(j-delta):(j+delta))); 
nuKmean(l  ,j)=mean(nuKfmal(l  ,(j-delta):(j+delta))); 

%calculate  the  atmospheric  density 
if  eKmagmean(lj)<0.1 

density  (j,p)--da(j,l).*1000./86400.*(BC).*(l./Vmean(lj).A2).*(mmKmean(jsl).*sqrt(l- 
eKmagmean(l  j)  A2)./sqrt(l  +eKmagmean(l  J).A2+2.*eKmagmean(l  j)  *cos(nuKmean(l  j)))); 
else 

density  (j,p)-0; 
end 
end 
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%plot 

figure; 

hold; 

scatter(Day(middle:count2-middle,l  ),density(middle:count2-middle,l )); 
xlabel(’time  [Day]’); 
ylabel('density  [kg/mA3]'); 

meandensity=mean(density(middle:count2- middle)); 

fprintf('\n  The  mean  density  using  %i  data  points  per  density  calculation  is  %\\  datapoints,  meandensity); 


